其他
一个完整的实证程序是什么样子, 以logit或ologit为例
凡是搞计量经济的，都关注这个号了
投稿：econometrics666@sina.cn
今天的主题：为了照固到经验不那么多的学者，以下是@甜菜 群友(大咖系列)分享给圈子的一个比较完整的实证程序, 以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 mx1`k'=mean(p`k') if x==1
egen mx0`k'=mean(p`k') if x==0
local k= `k'+1
}
**把预测概率值提出到scalar
forvalues l=0/1 {
forvalues j= 1/3 {
sum mx`l'`j'
scalar x`l'`j'=r(mean)
}
}
**把预测概率值显示出来
forvalues l=0/1 {
forvalues j= 1/3 {
dis as text "先生(女士), x`l'`j'=" as result x`l'`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 b`i'=el(ape,1,`i')
return scalar ape`i' =b`i'
}
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等相关模型
2.5年，计量经济圈近1000篇不重类计量文章，
可直接在公众号菜单栏搜索任何计量相关问题,
Econometrics Circle
数据系列：空间矩阵 | 工企数据 | PM2.5 | 市场化指数 | CO2数据 | 夜间灯光 | 官员方言 | 微观数据 | 内部数据计量系列：匹配方法 | 内生性 | 工具变量 | DID | 面板数据 | 常用TOOL | 中介调节 | 时间序列 | RDD断点 | 合成控制 | 200篇合辑 | 因果识别 | 社会网络 | 空间DID数据处理：Stata | R | Python | 缺失值 | CHIP/ CHNS/CHARLS/CFPS/CGSS等 |干货系列：能源环境 | 效率研究 | 空间计量 | 国际经贸 | 计量软件 | 商科研究 | 机器学习 | SSCI | CSSCI | SSCI查询 | 名家经验计量经济圈组织了一个计量社群，有如下特征：热情互助最多、前沿趋势最多、社科资料最多、社科数据最多、科研牛人最多、海外名校最多。因此，建议积极进取和有强烈研习激情的中青年学者到社群交流探讨，始终坚信优秀是通过感染优秀而互相成就彼此的。