ologit模型与logit_一个完整的实证程序, 以logit或ologit为例
// Copyrights belong to @计量经济圈
**For Econometrics Circle, 以下只是例子而已
use econometrics.dta, clear //输入数据
// Data management
global m F3A F4A A1 C18AA C18AB C18AC C18AD C18AE C18AF E8A B3A F6A //定义全局宏
**一个简单的循环来把变量的label全部删除
foreach varlist of global m {
label variable `varlist' ""
}
**数据整理基本要素
qui {
recode F3A (1=0)(3=1)(2=2)(else=.), gen(z)
label var z "wich one has higher z"
label define z 0 "a" 1 "b" 2 "c"
label values z
recode F4A (0=0)(1=1)(2=2)(3=3)(4=4)(5=5)(6=6)(7=7)(8=8)(10=10)(else=.), gen(number_child)
label var number_child "number of children"
recode A1 (1=0)(2=1)(else=.), gen(gender)
label var gender "if female"
gen Income = C18AA + C18AB + C18AC + C18AD + C18AE + C18AF
gen Inc = Income if inrange(Income,0,1000000)
label var Inc "total income"
capture drop Lninc
gen Lninc = ln(Inc+1)
gen income=Lninc
recode E8A (3=1) (1=0) (else=.), gen(x)
label var x "if x"
gen age2 = age^2
recode B3A (1=1) (2=2) (3=3) (4/5=4) (6=5) (7/8=6) (else=.), gen(education)
label var education"Education level"
recode F6A (1=0)(3=1)(2=2)(else=.), gen(y)
label var y "which one has higher y"
label define y 0 "a" 1 "b" 2 "c"
label values y y
}
/***********************************************************
Table 1: ordinal logit and generalized ordinal logit model
************************************************************/
drop if number_child<1 | number_child>6
global xlist x number_child education age age2 income z
global margef cttop(margeff) dec(3) alpha(0.01, 0.05, 0.1)e(ll r2_a r2_p) see replace
ologit y $xlist
estat ic // get AIC BIC
est store m
forval i = 0/2 {
est res m
margins, dydx(*) predict(outcome(`i')) post
est store m`i'
} //export marginal effect
esttab m0 m1 m2 using ologit1.rtf,pr2 b(%9.3f) star(* 0.1 ** 0.05 *** 0.01) se ar2 replace
//自己安装了gologit2之后把前面的星号去掉即可运行
**version 14.0
**gologit2 y $xlist, autofit
**margeff, at(mean) replace
**outreg2 using gologit.doc, $margef
**estat ic // get AIC BIC
/******************************************
Table 2: average predicted probability
*******************************************/
ologit y $xlist
local j=1
forvalues j=1(1)3 {
capture drop p`j'
capture drop mx1`j'
capture drop mx0`j'
local j=`j'+1
}
local k=1
forvalues k= 1(1)3{
predict p`k'
egen mx1k'=mean(pk') if x==1
egen mx0k'=mean(pk') if x==0
local k= `k'+1
}
**把预测概率值提出到scalar
forvalues l=0/1 {
forvalues j= 1/3 {
sum mxl'j'
scalar xl'j'=r(mean)
}
}
**把预测概率值显示出来
forvalues l=0/1 {
forvalues j= 1/3 {
dis as text "先生(女士), xl'j'=" as result xl'j'
}
}
/*************************************
Table 3: regression by subgroups
**************************************/
**按照教育和孩子数量分别回归
global sub "education number_child"
local k=1
qui {
foreach subgroup of global sub {
forvalues k=1/6{
ologit y $xlist if subgroup'==k'
forval i = 0/2 {
est res m
margins, dydx(*) predict(outcome(`i')) post
est store m`i'
}
outreg2 using education`k'.doc
}
}
}
**与前面的回归结果是一模一样的
by education, sort: ologit y $xlist
by number_child, sort: ologit y $xlist
/*************************************
Table 4: Two types of standard errors
**************************************/
**自助标准误
bootstrap,reps(50): ologit y $xlist, vce(robust)
**Delta方法标准误
capture prog drop brape
program brape,rclass
syntax varlist
probit `varlist'
margins, dydx(_all)
matrix ape=r(b)
local k: word count `varlist'
local k=`k'-1
forvalues i=1/`k' {
scalar bi'=el(ape,1,i')
return scalar apei' =bi'
}
end
brape y $xlist
/*************************************
Table 5: Three types of models for binary response models
**************************************/
- probit, logit and cloglog estimation and output
recode y (0/1=0)(2=1), gen(y1)
qui {
probit y1 $xlist
est store probit
margins, dydx(_all) post
est store probit1
logit y1 $xlist
est store logit
margins, dydx(_all) post // post is for saving margins
est store logit1
cloglog y1 $xlist
est store cloglog
margins, dydx(_all) post
est store cloglog1
}
est table probit logit cloglog,star(0.1 0.05 0.01) stats(N p_r2 aic bic) //系数比较
est table probit1 logit1 cloglog1,star(0.1 0.05 0.01) //边际效应比较
/*************************************
Table 6: Additional content for review
**************************************/
sysuse auto, clear
generate price2 = price[_n-1] //price2变量第二行的数等于price变量第一行的数
regress price mpg foreign rep78
display _b[mpg] //显示mpg回归的系数
**multinomial logit model
mlogit rep78 gear_ratio displacement foreign
display [5]_b[foreign]
webuse dollhill2, clear
by age (smokes), sort: generate wgt=pyears[_N] //按照age, smokes的顺序排序
bysort age: generate wgt=pyears[_N] //组内age生成wgt, 在组内都是这组最后一个数值
by age, sort: egen mean_w = mean(smokes) //能够按照age生成smokes的均值
**异方差probit估计的优越性hetprobit
webuse hetprobxmpl, clear
probit y x
predict xhet1, pr
gen residual=y-xhet1
line residual x //看看是不是同方差，如果有回归趋势那证明是异方差
hetprobit y x, het(xhet) //异方差probit估计
line xhet x //看看是不是同方差
**关注一下ivprobit, meologit, xtprobit, clogit, cloglog等相关模型