一个完整的实证程序是什么样子, 以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等相关模型