glm.ado

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

	if `"`fam'"'~="glim_v1" & `"`fam'"'~="glim_v2" & /*
	*/ `"`fam'"'~="glim_v3" & `"`fam'"'~="glim_v4" & /*
	*/ `"`fam'"'~="glim_v5" & `"`fam'"'~="glim_v6" {
		tokenize `"`f'"'
		if `"`2'"' != "" {
			global SGLM_fa `2'
		}
		if `"`3'"'!="" {
			di as err "family(`f') invalid"
			exit 198
		}
	}
	MapLink `ulink'
        local link `"`r(link)'"'
        local pow `"`r(power)'"'
        local scale1 = `"`fam'"'=="glim_v2" | /*
			*/ `"`fam'"'=="glim_v3" | `"`fam'"'=="glim_v6"

        if `"`link'"'=="" {             /* apply defaults */
                local link "glim_l11"
                if `"`fam'"'=="glim_v1"      { local link "glim_l01"
					       local pow 1           }
                else if `"`fam'"'=="glim_v2" { local link "glim_l02"
					       local pow 0           }
                else if `"`fam'"'=="glim_v3" { local link "glim_l03"
					       local pow  0          }
                else if `"`fam'"'=="glim_v4" { local link "glim_l09"
					       local pow -1          }
                else if `"`fam'"'=="glim_v5" { local link "glim_l10"
					       local pow -2          }
                else if `"`fam'"'=="glim_v6" { local link "glim_l03"
					       local pow  0          }
                else local link
        }

        ret local family `"`fam'"'
        ret local link   `"`link'"'
        ret local power  `pow'
        ret local scale  `scale1'
        ret local m      `m'
        ret local mfixed `mfixed'
        ret local k      `k'
end


program MapFam, rclass /* <code> */
	args f

        local s1
        local f = lower(trim(`"`f'"'))
        local l = length(`"`f'"')

        if `"`f'"'==""                                    { local s1 "1" }
        else if `"`f'"'==substr("gaussian",1,max(`l',3))  { local s1 "1" }
        else if `"`f'"'==substr("normal",1,`l')           { local s1 "1" }
        else if `"`f'"'==substr("binomial",1,`l')         { local s1 "2" }
        else if `"`f'"'==substr("bernoulli",1,`l')        { local s1 "2" }
        else if `"`f'"'==substr("poisson",1,`l')          { local s1 "3" }
        else if `"`f'"'==substr("gamma",1,max(`l',3))     { local s1 "4" }
        else if `"`f'"'==substr("igaussian",1,max(`l',2)) { local s1 "5" }
        else if `"`f'"'==substr("inormal",1,max(`l',2))   { local s1 "5" }
        else if `"`f'"'=="ivg"                            { local s1 "5" }
        else if `"`f'"'==substr("nbinomial",1,max(2,`l')) { local s1 "6" }

	if "`s1'" != "" {
		local s1 "glim_v`s1'"
	}
	else {
		local s1 "`f'"
	}
        ret local famcode `s1'
end

program MapLink, rclass /* <code> <coderest> ... */
	args ulink upow rest
        if `"`rest'"'!="" {
                di as err "link(`ulink' `upow' `rest') invalid"
                exit 198
        }

        local ulink = lower(`"`ulink'"')
        local l     = length(`"`ulink'"')
        local s1      "11"                /* power(a) link */
        local s2      .                   /* a of power(a) */

        if `"`ulink'"'==""                              { local s1 "" }
        else if `"`ulink'"'==substr("identity",1,`l')   { local s1 "01"
							  local s2 1    }
        else if `"`ulink'"'==substr("reciprocal",1,`l') { local s1 "09"
							  local s2 -1   }
        else if `"`ulink'"'=="log"                      { local s1 "03"
							  local s2 0    }
        else if `"`ulink'"'==substr("power",1,max(`l',3)) {
                capture confirm number `upow'
                if _rc {
                        di as err "invalid # in link(power #)"
                        exit 198
                }
		if `upow' == 0 {
			local s1 "03"
		}
		if `upow' == -1 {
			local s1 "09"
		}
		else if `upow' == -2 {
			local s1 "10"
		}
                local s2 `upow'
                local upow
        }
        else if `"`ulink'"'==substr("opower",1,max(`l',3)) {
                capture confirm number `upow'
                if _rc {
                        di as err "invalid # in link(opower #)"
                        exit 198
                }
                if `upow'==0 {
                        local s1 "02"	/* logit */
                }
                else {
                        local s1 "12"   /* odds-power */
                }
		local s2 `upow'
                local upow
        }
        else if `"`ulink'"'==substr("logit",1,`l')          { local s1 "02"
							      local s2 0    }
        else if `"`ulink'"'==substr("nbinomial",1,max(`l',2)) { local s1 "04" }
        else if `"`ulink'"'==substr("logc",1,4)             { local s1 "05" }
        else if `"`ulink'"'==substr("loglog",1,max(`l',4))  { local s1 "06" }
        else if `"`ulink'"'==substr("cloglog",1,`l')        { local s1 "07" }
        else if `"`ulink'"'==substr("probit",1,`l')         { local s1 "08" }

        else {
		local s1
		local s2 "`upow'"
        }

	if `"`s1'"' != "" {
		local s1 "glim_l`s1'"
	}
	else {
		local s1 "`ulink'"
	}

        ret local link  `s1'
        ret local power `s2'
end

program Eform , rclass
        local var "`e(varfunc)'"
        local lnk "`e(link)'"

        local eform "ExpB"
        if ("`var'" == "glim_v3" | "`var'" == "glim_v6") {
                if "`lnk'" == "glim_l03" {
                        local eform "IRR"
                }
        }
        else if "`var'" == "glim_v2" {
                if "`lnk'" == "glim_l02" {
                        local eform "Odds Ratio"
                }
                else if "`lnk'" == "glim_l03" {
                        local eform "Risk Ratio"
                }
                else if "`lnk'" == "glim_l05" {
                        local eform "Hlth Ratio"
                }
                else local eform "ExpB"
        }
        return local eform "`eform'"
end

program HAC, rclass
	args hacnam haclag wv res touse constant xvars

	local xv "`xvars'"
	if "`constant'" == "" {
		tempvar CONS
		gen byte `CONS' = 1
		local xn "`xv' _cons"
		local xv "`xv' `CONS'"
	}
	else {  local xn "`xv'" }
	local nx : word count `xv'

	qui summ `touse', meanonly
	local nobs = r(N)

	tempvar vt1
	tempname tt xtx tx

	capture tsset
	if _rc == 0 {
		sort `r(panelvar)' `r(timevar)'
	}

	gen double `vt1' = .
	local j 0
	while `j' <= `haclag' {
		capture mat drop `tt'
		local i 1
		while `i' <= `nx' {
			local x : word `i' of `xv'
			tsrevar `x'
			local x `r(varlist)'
			qui replace `vt1' = `x'[_n-`j']*`res'*`res'[_n-`j']* /*
					*/ `wv'[_n-`j'] if `touse'
			mat vecaccum `tx' = `vt1' `xv' if `touse', nocons
			mat `tt' = nullmat(`tt') \ `tx'
			local i = `i'+1
		}
		`hacnam' `haclag' `j'
		if `r(wt)'>=. {
			noi di as err "Newey-West weights are missing"
			exit 199
		}
		mat `tt' = (`tt' + `tt'')*r(wt)
		if `j' > 0 {
			mat `xtx' = `xtx'+`tt'
		}
		else {
			mat `xtx' = `tt'*.5
		}
		local j = `j'+1
	}
	mat colnames `xtx' = `xn'
	mat rownames `xtx' = `xn'
	return matrix V `xtx'
end

program LoadXi
	args beta xi i

	local p = colsof(`beta')
	local xv : colnames(`beta')
	local j 1
	while `j'<=`p' {
		local x : word `j' of `xv'
		if "`x'" != "_cons" {
			mat `xi'[1,`j'] = `x'[`i']
		}
		else    mat `xi'[1,`j'] = 1
		local j = `j'+1
	}
end

program JK1nr
	args b V y xvars eta mu v dmu tvar z moffset constant wt touse weight
	capture drop `eta'
	capture drop `mu'
	capture drop `v'
	capture drop `dmu'
	_predict `eta', xb
	$SGLM_L 1 `eta' `mu'
	$SGLM_V 1 `eta' `mu' `v'
	$SGLM_L 2 `eta' `mu' `dmu'
	capture drop `z'
	tempvar W h
	tempname ll
	gen double `z' = `eta' + (`y'-`mu')/`dmu' `moffset' if `touse'
	gen double `W' = `dmu'*`dmu'/`v' if `touse'
	estimates hold `ll'
	cap noisily quietly regress `z' `xvars' [iw=`W'*`wt'], mse1 `constant'
	if e(N)>=. {exit 459}
	if _rc {exit _rc}
	_predict double `h' if `touse', stdp
	replace `h' = `h'*`h'*`v'
	estimates unhold `ll'

	capture drop `z'
	gen double `z' = (`mu'-`y')/(1-`h')
	sort `touse' `tvar'
	capture drop `dmu'

	gen `dmu' = _n if `touse'
	local p = colsof(`b')

	tempname bi xtxi xi ss
	mat `xi' = J(1,`p',0)
	mat `V' = 0*`V'
	if "`weight'" == "fweight" {
		mat accum `xtxi' = `xvars' [iw=`v'*`wt'] if `touse', `constant'
	}
	else    mat accum `xtxi' = `xvars' [iw=`v'] if `touse', `constant'
	mat `xtxi' = syminv(`xtxi')
	summ `dmu'
	local factor 1
	local max = r(max)
	local i   = r(min)
	while `i' <= `max' {
		if !`touse'[`i'] {
			local i = `i' + 1
			continue
		}
		LoadXi "`b'" "`xi'" `i'
		mat `bi' = `xi'*`xtxi'
		scalar `ss' = `z'[`i']
		mat `bi' = `ss'*`bi'
		if "`weight'" == "fweight" { local factor = `wt'[`i'] }
		mat `V' = `V' + `factor'*(`bi')'*`bi'
		local i = `i'+1
	}
end

program JK1irls
	args b V y xvars h eta mu v dmu tvar z moffset constant wt touse weight
	capture drop `v'
	capture drop `dmu'
	capture drop `z'
	$SGLM_V 1 `eta' `mu' `v'
	gen double `z' = (`mu'-`y')/(1-`h')
	sort `touse' `tvar'
	gen `dmu' = _n if `touse'
	local p = colsof(`b')

	tempname bi xtxi xi ss
	mat `xi' = J(1,`p',0)
	mat `V' = 0*`V'
	mat accum `xtxi' = `xvars' [iw=`v'*`wt'] if `touse', `constant'
	mat `xtxi' = syminv(`xtxi')
	summ `dmu'
	local factor 1
	local max = r(max)
	local i = r(min)
	while `i' <= `max' {
		if !`touse'[`i'] {
			local i = `i' + 1
			continue
		}
		LoadXi "`b'" "`xi'" `i'
		mat `bi' = `xi'*`xtxi'
		scalar `ss' = `z'[`i']
		mat `bi' = `ss'*`bi'
		if "`weight'" == "fweight" { local factor = `wt'[`i'] }
		mat `V' = `V' + `factor'*(`bi')'*`bi'
		local i = `i'+1
	}
end

program JKnife
	args b V y xvars tvar offset constant family irls link /*
		*/ wtopt touse cluster dots weight wt

	if "`dots'" != "" { noi di }
	sort `touse' `tvar'
	local p = colsof(`b')
	mat `V' = 0*`V'

	if "`irls'" == "" {
		local iopt "from(`b')"
	}
	else {
		tempvar xb mu
		qui _predict double `xb' if `touse', xb
		qui $SGLM_L 1 `xb' `mu'
		local iopt "mu(`mu')"
	}
	local reg = "$SGLM_V" == "glim_v1" & "$SGLM_L" == "glim_l01" & "`irls'" == ""
	if `reg' { local iopt }
	tempname bi ll
	estimates hold `ll'
	tempvar dd
	if "`cluster'" != "" {
		qui egen `dd' = group(`cluster') if `touse'
		qui summ `dd'
		local nsize = r(max)
	}
	else    {
		qui gen  `dd' = _n if `touse'
		qui summ `dd'
		local nsize = r(N)
	}
	local max = r(max)
	local min = r(min)
	local i   = `min'
	local nb = 0

	if "`dots'" != "" {
		di as txt "Jackknife iterations (" /*
			*/ as res `nsize' as txt ")"
		di as txt "{hline 4}{c +}{hline 3} 1 "/*
			*/"{hline 3}{c +}{hline 3} 2 "/*
			*/"{hline 3}{c +}{hline 3} 3 "/*
			*/"{hline 3}{c +}{hline 3} 4 "/*
			*/"{hline 3}{c +}{hline 3} 5"
	}

	if "`weight'" == "fweight" {
		local wtopt "[fweight=`wt']"
		local fixwt 1
	}
	else    local fixwt 0


	local nx : word count `xvars'

	while `i' <= `max' {
		if "`cluster'" == "" & !`touse'[`i'] {
			local i = `i'+1
			continue
		}
		if `fixwt' == 0 {
			qui _rmcoll `xvars' `wtopt' if `touse' & `dd'!=`i', `constant'
			local nxrm : word count `r(varlist)'
		}
		else    local nxrm = `nx'

		if `nxrm' == `nx' {
			capture {
				if `fixwt' {
					qui summ `wt' if `dd'==`i', meanonly
					local factor = r(mean)
					qui replace `wt'=`wt'-1 if `dd'==`i'
					qui glm_7 `y' `xvars' if `touse' /*
						*/ `wtopt' /*
						*/ , `family' /*
						*/ `link' `offset' /*
						*/ `constant' `irls' `iopt'
				}
				else {
					local factor 1
					qui glm_7 `y' `xvars' if `touse' & /*
					*/ `dd'!=`i' `wtopt' /*
					*/ , `family' /*
					*/ `link' `offset' /*
					*/ `constant' `irls' `iopt'
				}
				mat `bi' = e(b)
				mat `V' = `V' + `factor'*(`bi'-`b')'*(`bi'-`b')
			}
			if _rc {
				if _rc == 1 { exit _rc }
				local err 1
				local nb = `nb'+`factor'
			}
			else    local err 0
			if `fixwt' {
				qui replace `wt'=`wt'+1 if `dd'==`i'
			}
		}
		else    local err 1
		if "`dots'" != "" {
			local j = `i'-`min'+1
			if `err' {
				noi di as err "x" _c
			}
			else {
				noi di as txt "." _c
			}
			if int(`j'/50)*50 == `j' {
				noi di as res " " `j'
			}
		}
		local i = `i'+1
	}
	estimates unhold `ll'
	global SGLM_nc = `nsize'
	global SGLM_bm = `nb'
	if "`dots'" != "" { noi di }
end

program Bstrap
	args b V y xvars offset constant family irls link /*
		*/ wtopt touse brep clopt dots

	if "`dots'" != "" { noi di }
	sort `touse' `tvar'
	local p = colsof(`b')
	mat `V' = 0*`V'

	if "`irls'" == "" {
		local iopt "from(`b')"
	}
	else {
		tempvar xb mu
		qui _predict double `xb' if `touse', xb
		qui $SGLM_L 1 `xb' `mu'
		local iopt "mu(`mu')"
	}

	local reg = "$SGLM_V" == "glim_v1" & "$SGLM_L" == "glim_l01" & "`irls'" == ""
	if `reg' { local iopt }

	tempname bi ll
	preserve
	qui drop if `touse'==0
	tempfile tmp
	estimates hold `ll'
	qui save `tmp'
	local nb = 0
	local i 1

	if "`dots'" != "" {
		di as txt "Bootstrap iterations (" /*
			*/ as res `brep' as txt ")"
		di as txt "{hline 4}{c +}{hline 3} 1 "/*
                        */"{hline 3}{c +}{hline 3} 2 "/*
                        */"{hline 3}{c +}{hline 3} 3 "/*
                        */"{hline 3}{c +}{hline 3} 4 "/*
                        */"{hline 3}{c +}{hline 3} 5"
	}

	local nx : word count `xvars'

	while `i' <= `brep' {
		qui use `tmp', clear
		qui bsample , `clopt'

		qui _rmcoll `xvars' `wtopt', `constant'
		local nxrm : word count `r(varlist)'

		local err 0
		if `nxrm' == `nx' {
			capture {
				qui glm_7 `y' `xvars' `wtopt' /*
					*/ , `family' /*
					*/ `link' `offset' /*
					*/ `constant' `irls' `iopt'
				mat `bi' = e(b)
				mat `V' = `V' + (`bi'-`b')'*(`bi'-`b')
			}
			if _rc {
				local err 1
				if _rc == 1 { exit _rc }
				local nb = `nb'+1
			}
			else    local err 0
		}
		else local err 1
		if "`dots'" != "" {
			if `err' {
				noi di as err "x" _c
			}
			else {
				noi di as txt "." _c
			}
			if int(`i'/50)*50 == `i' {
				noi di as res " " `i'
			}
		}
		local i = `i'+1
	}
	restore
	estimates unhold `ll'
	global SGLM_bm = `nb'
	if "`dots'" != "" { noi di }
end
exit