一个完整的实证程序, 以logit或ologit为例

凡是搞计量经济的，都关注这个号了

投稿：[e-mail]

所有计量经济圈方法论丛的程序文件, 微观数据库和各种软件都放在社群里.欢迎到计量经济圈社群交流访问.

题外话：有些圈友本身基础比较差，以为报个培训班就能够掌握计量方法以及数据管理软件，其实是比较危险的想法。在计量经济圈社群里的大咖，谁不是经过各种问题的洗礼而增进一点点理解的？如果以下这些Tips和Tricks都掌握了，你已经属于高手中的高手了。

10.为发表SSCI刊物, 必看这21篇文章(人大复印资料转载)

今天的主题：为了照顾到经验不那么多的学者，以下是@甜菜 群友(大咖系列)分享给圈子一个比较完整的实证程序, 以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等相关模型

为了照顾到群外的圈友，因此转发5分钟后截图到公众号后台可以得到里面的数据和程序do文档。

2年，计量经济圈公众号近1000篇文章，

Econometrics Circle

数据系列：空间矩阵 | 工企数据 | PM2.5 | 市场化指数 | CO2数据 |  夜间灯光 

计量系列：匹配方法 | 内生性问题 | 工具变量 | DID |  面板数据 | 常用TOOL

计量经济圈组织了一个计量社群，有如下特征：热情互助最多、前沿趋势最多、社科资料最多、社科数据最多、科研牛人最多、海外名校最多。因此，建议积极进取和有强烈研习激情的中青年学者到社群交流探讨，始终坚信优秀是通过感染优秀而互相成就彼此的。