xtgee.ado

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

			mat `wt'[`j',`i'] = `tmp'[1,`k']
			mat `wt'[`i',`j'] = `tmp'[1,`k']
			local j = `j'+1
		}
		local i = `i'+1
	}

	mat `wt' = syminv(`wt')
	mat `al' = J(1,`b2',0)
	local i 1
	while `i' <= `b2' {
		local j = `i'+1
		mat `al'[1,`i'] = `a'[1,`j']
		local i = `i'+1
	}
	mat `al' = `al'*`wt'
	mat drop `wt'

	tempname scnm
	local i = `b'
	while `i' <= `n' {
		scalar `scnm' = 0
		local j 1
		while `j' <= `b2' {
			local k = `i'-`j'
			scalar `scnm' = `scnm' + `al'[1,`j']*`tmp'[1,`k']
			local j = `j'+1
		}
		mat `tmp'[1,`i'] = `scnm'
		local i = `i'+1
	}
	_GEERC , c(arm) r(`R') n(`n') a(`tmp') b(`b')
end


program define checkcm
	args cmat

	cap local nrow = rowsof(`cmat')
	if _rc {
		noi di in red "`cmat' not a matrix"
		exit 198
	}

	local ncol = colsof(`cmat')

	if `ncol' != `nrow' {
		noi di in red "correlation matrix not square"
		exit 198
	}

	local i 1
	while `i' < `nrow' {
		if `cmat'[`i',`i'] != 1.0 {
			noi di in red "`cmat' not a correlation matrix"
			exit 198
		}
		local j = `i' + 1
		while `j' <= `nrow' {
			if `cmat'[`i',`j'] != `cmat'[`j',`i'] {
				noi di in red "`cmat' not symmetric"
				exit 198
			}
			if `cmat'[`i',`j'] > 1 | `cmat'[`i',`j'] < -1 {
				noi di in red "`cmat' not a correlation matrix"
				exit 198
			}
			local j = `j'+1
		}
		local i = `i'+1
	}
	if `cmat'[`i',`i'] != 1.0 {
		noi di in red "`cmat' not a correlation matrix"
		exit 198
	}
end


program define checkfam /* family link */
	args family link

	local flag 0

	if "`link'" == "power" {
		local flag 0
	}

	else if "`family'"=="binom" {
		if "`link'"=="nbinom" {
			local flag 1
		}
	}
	else if "`family'"=="gauss" {
		if "`link'"!="ident" & "`link'"!="log"  & "`link'"!="recip" {
			local flag 1
		}
	}

	else if "`family'"=="poiss" {
		if "`link'"!="ident"&"`link'"!="log" & "`link'"!="recip" {
			local flag 1
		}
	}

	else if "`family'" == "gamma" {
		if "`link'"!="ident" & "`link'"!="log" & "`link'" != "recip" {
			local flag 1
		}
	}

	else if "`family'" == "nbinom" {
		if "`link'"!="ident" & "`link'"!="log" & "`link'"!="nbinom"{
			local flag 1
		}
	}

	else if "`family'" == "igauss" {
		if "`link'"!="ident" & "`link'"!="log" {
			local flag 1
		}
	}

	if `flag' == 1 {
		noi di in red "unsupported family-link combination"
		exit 198
	}
end

program define GetName
	args link fam corr

	global S_1	/* link		*/
	global S_2	/* family	*/
	global S_3	/* corr		*/
	global S_4	/* family arg	*/
	global S_5 -1	/* corr lag	*/
	global S_6	/* corr mat	*/

	local l = length("`link'")
	local f = lower(trim("`link'"))
	local f : word 1 of `f'
	local l = length("`f'")

	if "`f'" == substr("identity",1,max(`l',1))     { global S_1 "ident" }
	else if "`f'" == substr("log",1,max(`l',3))     { global S_1 "log"   }
	else if "`f'" == substr("logit",1,max(`l',4))   { global S_1 "logit" }
	else if "`f'" == substr("power",1,max(`l',3))   { global S_1 "power" }
	else if "`f'" == substr("opower",1,max(`l',3))  { global S_1 "opower" }
	else if "`f'" == substr("cloglog",1,max(`l',2)) { global S_1 "cloglog" }
	else if "`f'" == substr("probit",1,max(`l',1))  { global S_1 "probit"}
	else if "`f'" == substr("nbinomial",1,max(`l',2)){ global S_1 "nbinom"}
	else if "`f'"==substr("reciprocal",1,max(`l',3)){ global S_1 "recip" }
	else if "`f'" != "" {
		noi di in red "unknown link() `link'"
		exit 198
	}

	global S_X_pow
	if "$S_1" == "power" | "$S_1" == "opower" {
		tokenize `link'
		global S_X_pow "`2'"
		if "$S_X_pow" == "" { global S_X_pow 1 }
		capture confirm number $S_X_pow
		if _rc {
			noi di in red "$S_1 link requires numeric argument"
			exit 198
		}
	}

	local l = length("`fam'")
	local f = lower(trim("`fam'"))
	local f : word 1 of `f'
	local l = length("`f'")

	if "`f'" == substr("gaussian",1,max(`l',3))      { global S_2 "gauss" }
	else if "`f'" == substr("normal",1,max(`l',3))   { global S_2 "gauss" }
	else if "`f'" == substr("igaussian",1,max(`l',2)){ global S_2 "igauss"}
	else if "`f'" == substr("poisson",1,max(`l',1))  { global S_2 "poiss" }
	else if "`f'" == substr("nbinomial",1,max(`l',2)){ global S_2 "nbinom"}
	else if "`f'" == substr("binomial",1,max(`l',1)) { global S_2 "binom" }
	else if "`f'" == substr("gamma",1,max(`l',3))    { global S_2 "gamma" }
	else if "`f'" != "" {
		noi di in red "unknown family() `fam'"
		exit 198
	}
	
	tokenize `fam'
	global S_4 "`2'"
	global S_X_nba
	if "$S_2" == "nbinom" {
		if "$S_4" == "" { global S_4 1 }
		capture confirm number $S_4
		if _rc{  
			noi di in red /*
			*/ "family(nbinomial #) requires # to be positive"
			exit 198
		}
		else if $S_4 <=0 {
			noi di in red /*
			*/ "family(nbinomial #) requires # to be positive"
			exit 198
		}
		global S_X_nba = 1/$S_4
		global S_4
	}

	GetCorr "`corr'" /* sets S_3, S_5, S_6 */

	if "$S_3" == "" {
		global S_3 "exchan"
	}

	if "$S_2" == "" {
		global S_2 "gauss"
	}

	if "$S_1" == "" {
		if "$S_2" == "gauss"  { global S_1 "ident" }
		if "$S_2" == "poiss"  { global S_1 "log"   }
		if "$S_2" == "binom"  { global S_1 "logit" }
		if "$S_2" == "gamma"  { global S_1 "recip" }
		if "$S_2" == "igauss" {
			global S_1 "power"
			global S_X_pow = -2
		}
		if "$S_2" == "nbinom" { global S_1 "log"   }
	}

end

program define GetCorr /* correlation <#> */
/*
	We parse
		<correlation>
		<correlation>[ ]#
		<correlation> <matname>
*/
	ReTok `*'
	local f = lower(trim("$S_3"))
	local l = length("`f'")
	local rest "$S_5"

	global S_3	/* corr */
	global S_5 -1	/* lag */
	global S_6	/* matname */

	if "`f'"=="" { exit }

	if "`f'" == substr("independent",1,max(`l',3))	     {
		global S_3 "indep"
		NoLag `f' `rest'
	}
	else if "`f'" == substr("arm",1,max(`l',2))   {
		global S_3 "arm"
		SetLag `f' `rest'
	}
	else if "`f'" == substr("exchangeable",1,max(`l',3))   {
		global S_3 "exchan"
		NoLag `f' `rest'
	}
	else if "`f'" == substr("statm",1,max(`l',3))   {
		global S_3 "statm"
		SetLag `f' `rest'
	}
	else if "`f'" == substr("stationary",1,max(`l',3))   {
		global S_3 "statm"
		SetLag `f' `rest'
	}
	else if "`f'" == substr("nonstationary",1,max(`l',3))   {
		global S_3 "nonst"
		SetLag `f' `rest'
	}
	else if "`f'" == substr("unstructured",1,max(`l',3))   {
		global S_3 "unstr"
		NoLag `f' `rest'
	}
	else if "`f'" == substr("fixed",1,max(`l',3))   {
		global S_3 "fixed"
		SetMat `f' `rest'
		checkcm $S_6
	}
	else if "`f'" != "" {
		noi di in red "unknown corr() structure `corr'"
		exit 198
	}
end

program define ReTok
	if "`2'"=="" {
		local c = substr("`1'",length("`1'"),1)
		while "`c'">="0" & "`c'"<="9" {
			local 1 = substr("`1'",1,length("`1'")-1)
			local 2 "`c'`2'"
			local c = substr("`1'",length("`1'"),1)
		}
	}
	global S_3 "`1'"
	mac shift
	global S_5 "`*'"
end

program define NoLag
	if "`2'"!="" {
		BadCorr `*'
	}
end

program define BadCorr
	di in red "corr(`*') invalid"
	exit 198
end

program define SetLag
	if "`3'" != "" {
		BadCorr `*'
	}
	if "`2'"=="" {
		global S_5 1
		exit
	}
	cap confirm integer number `2'
	if _rc {
		BadCorr `*'
	}
	if `2'<1 {
		BadCorr `*'
	}
	global S_5 `2'
end

program define SetMat
	if "`3'"!="" {
		BadCorr `*'
	}
	cap mat list `2'
	if _rc {
		if _rc==111 {
			di in red "matrix `3' not found"
			exit 111
		}
		BadCorr `*'
	}
	global S_6 `2'
end

program define FixScale
	args y mu m wtexp

	global S_X_fix
	if "$S_E_scp" == "" {
		global S_X_fix 1
		if "$S_X_mvar" == "gauss" | "$S_X_mvar" == "gamma" | /*
			*/ "$S_X_mvar" == "igauss" {
			global S_E_scp "x2"
		}
		else {
			global S_E_scp 1
		}
	}

	tempvar dj
	if "$S_X_mvar" == "gauss" {
		if "$S_E_scp" == "dev" {
			gen double `dj' = (`y'-`mu')^2
			global S_X_fix 1
		}
		else if "$S_E_scp" == "x2" {
			gen double `dj' = (`y'-`mu')
		}
	}
	else if "$S_X_mvar" == "igauss" {
		if "$S_E_scp" == "dev" {
			gen double `dj' = (`y'-`mu')^2 / (`mu'^2*`y')
			global S_X_fix 1
		}
		else if "$S_E_scp" == "x2" {
			gen double `dj' = (`y'-`mu') / sqrt(`mu'^3)
		}
	}
	else if "$S_X_mvar" == "gamma" {
		if "$S_E_scp" == "dev" {
			gen double `dj' = -2*(log(`y'/`mu')-(`y'-`mu')/`mu')
			global S_X_fix 1
		}
		else if "$S_E_scp" == "x2" {
			gen double `dj' = (`y'-`mu')/sqrt((`mu')^2)
		}
	}
	else if "$S_X_mvar" == "poiss" {
		if "$S_E_scp" == "dev" {
			gen double `dj' = 2*`mu'
			replace `dj' = 2*(`y'*log(`y'/`mu')-(`y'-`mu')) if `y'
			global S_X_fix 1
		}
		else if "$S_E_scp" == "x2" {
			gen double `dj' = (`y'-`mu')/sqrt(`mu')
			global S_X_fix 1
		}
	}
	else if "$S_X_mvar" == "binom" {
		if "$S_E_scp" == "dev" {
			gen double `dj' = 2*`y'*log(`y'/`mu')+ /*
			*/ 2*(`m'-`y')*log((`m'-`y')/(`m'-`mu'))
			replace `dj' = 2*`m'*log(`m'/(`m'-`mu')) if `y'==0
			replace `dj' = 2*`y'*log(`y'/`mu') if `y'==`m'
			global S_X_fix 1
		}
		else if "$S_E_scp" == "x2" {
			gen double `dj' = (`y'-`mu')/sqrt(`mu'*(1-`mu'/`m'))
			global S_X_fix 1
		}
	}
	else if "$S_X_mvar" == "nbinom" {
		local k = $S_X_nba
		if "$S_E_scp" == "dev" {
			gen double `dj' = 2*ln(1+`k'*`mu')/`k'
			replace `dj' = 2*(`y'*ln(`y'/`mu')-(1+`k'*`y')/`k'* /*
				*/ ln((1+`k'*`y')/(1+`k'*`mu'))) if `y'
			global S_X_fix 1
		}
		else if "$S_E_scp" == "x2" {
			gen double `dj' = (`y'-`mu')/sqrt(`mu'+`k'*`mu'^2)
			global S_X_fix 1
		}
	}

	if "$S_E_scp" == "x2" {
		replace `dj' = (`dj')^2
	}

	if "$S_E_scp" != "phi" {
		cap confirm number $S_E_scp
		if _rc {
			summ `dj' `wtexp'
			global S_1 = r(sum)/$S_E_dof
			global S_E_dis = $S_1
		}
		else {
			global S_1 = $S_E_scp
		}
		global S_X_fix 1
	}
        else    global S_1


	if "$S_X_mvar" == "gauss" {
		cap drop `dj'
		gen double `dj' = (`y'-`mu')^2
		summ `dj' `wtexp', meanonly
		global S_X_D1 = r(sum)
		global S_X_D2 = r(sum)/$S_E_dof
		cap drop `dj'
		gen double `dj' = (`y'-`mu')
		replace `dj' = (`dj')^2
		summ `dj' `wtexp'
		global S_X_P1 = r(sum)
		global S_X_P2 = r(sum)/$S_E_dof
	}
	else if "$S_X_mvar" == "igauss" {
		cap drop `dj'
		gen double `dj' = (`y'-`mu')^2 / (`mu'^2*`y')
		summ `dj' `wtexp', meanonly
		global S_X_D1 = r(sum)
		global S_X_D2 = r(sum)/$S_E_dof
		cap drop `dj'
		gen double `dj' = (`y'-`mu') / sqrt(`mu'^3)
		replace `dj' = (`dj')^2
		summ `dj' `wtexp', meanonly
		global S_X_P1 = r(sum)
		global S_X_P2 = r(sum)/$S_E_dof
	}
	else if "$S_X_mvar" == "poiss" {
		cap drop `dj'
		gen double `dj' = 2*`mu'
		replace `dj' = 2*(`y'*log(`y'/`mu')-(`y'-`mu')) if `y'
		summ `dj' `wtexp', meanonly
		global S_X_D1 = r(sum)
		global S_X_D2 = r(sum)/$S_E_dof
		cap drop `dj'
		gen double `dj' = (`y'-`mu')/sqrt(`mu')
		replace `dj' = (`dj')^2
		summ `dj' `wtexp', meanonly
		global S_X_P1 = r(sum)
		global S_X_P2 = r(sum)/$S_E_dof
	}
	if "$S_X_mvar" == "gamma" {
		cap drop `dj'
		gen double `dj' = -2*(log(`y'/`mu')-(`y'-`mu')/`mu')
		summ `dj' `wtexp', meanonly
		global S_X_D1 = r(sum)