. sysuse "nlsw88.dta", clear
. gen agesq = age*age
*-分组虚拟变量
. drop if race==3
. gen black = 2.race
. tab black
*-删除缺漏值
. global xx "ttl_exp married south hours tenure age* i.industry"
. reg wage $xx i.race
. keep if e(sample)
*-分组回归
. global xx "ttl_exp married south hours tenure age* i.industry"
. reg wage $xx if black==0
. est store White
. reg wage $xx if black==1
. est store Black
*-结果对比
. local m "White Black"
. esttab `m', mtitle(`m') b(%6.3f) nogap drop(*.industry) ///
s(N r2_a) star(* 0.1 ** 0.05 *** 0.01)
结果:
Table1: 白人组和黑人组工资影响因素差异对比
------------Table 1-------------------
(1) (2)
White Black
--------------------------------------
ttl_exp 0.251*** 0.269***
(6.47) (4.77)
married -0.737** 0.091
(-2.31) (0.23)
... (output ommited) ...
--------------------------------------
N 1615.000 572.000
r2_a 0.112 0.165
--------------------------------------
t statistics in parentheses
* p<0.1, ** p<0.05, *** p<0.01
dropvars ttl_x_black marr_x_black
global xx "ttl_exp married south hours tenure age* i.industry" //Controls
gen ttl_x_black = ttl_exp*black //交乘项
reg wage black ttl_x_black $xx //全样本回归+交乘项
为节省篇幅,仅列出最关键的结果如下:
. reg wage black ttl_x_black $xx //全样本回归+交乘
Number of obs = 2,187
R-squared = 0.1379
Adj R-squared = 0.1299
--------------------------------------------------
wage | Coef. SE t P>|t|
-------------+------------------------------------
black | -0.819 0.802 -1.02 0.307
ttl_x_black | -0.018 0.059 -0.31 0.756
ttl_exp | 0.254 0.035 7.23 0.000
married | -0.465 0.253 -1.84 0.066
... (output omitted) ...
--------------------------------------------------
*-Step1: 分别针对两个样本组执行估计
reg wage $xx if black==0
est store w //white
reg wage $xx if black==1
est store b //black
*-Step 2: SUR
suest w b
*-Step 3: 检验系数差异
test [w_mean]ttl_exp = [b_mean]ttl_exp
test [w_mean]married = [b_mean]married
test [w_mean]south = [b_mean]south
Step 2 的结果如下(为便于阅读,部分变量的系数未呈现):
. suest w b
Simultaneous results for w, b
Number of obs = 2,187
---------------------------------------------------
| Robus
| Coef. Std. Err. z P>|z|
------------+--------------------------------------
w_mean
ttl_exp | 0.251 0.036 6.92 0.000
married | -0.737 0.349 -2.11 0.035
south | -0.813 0.289 -2.81 0.005
... (output omitted) ..
------------+--------------------------------------
b_mean
ttl_exp | 0.269 0.051 5.25 0.000
married | 0.091 0.405 0.23 0.822
south | -2.041 0.425 -4.80 0.000
... (output omitted) ..
------------------+--------------------------------
对上述命令和结果的简要解释如下:
白人组和黑人组的估计结果分别存储于 w 和 b 两个临时性文件中;
执行 suest w b 命令时,白人组和黑人组的被视为两个方程,即文的 (2a) 和 (2b) 式。Stata 会自动将两个方程对应的样本联合起来,采用 GLS 执行似无相关估计(SUR);
使用 suest 时,允许两个样本组的解释变量个数不同。但由于一些技术上的问题尚未解决(很快可以解决掉),bdiff 命令要求两个样本组中的解释变量个数相同。在上例中,白人组在 Mining 行业的观察值个数为零(输入tab industry black 可以查看),导致我们加入行业虚拟变量时,白人组只有 10 个行业虚拟变量,而黑人组则有 11 个行业虚拟变量。为此,在上述命令中,我使用 drop if industry==2 命令删除了 Mining 行业的观察值。
*-SUEST test for panel data
*-数据概况
webuse "nlswork", clear
xtset idcode year
xtdes
*-对核心变量执行组内去心:去除个体效应
help center //外部命令, 下载命令为 ssc install center, replace
local y "ln_wage"
local x "hours tenure ttl_exp south"
bysort id: center `y', prefix(cy_) //组内去心
bysort id: center `x', prefix(cx_)
*-分组回归分析
reg cy_* cx_* i.year if collgrad==0 // 非大学组
est store Yes
reg cy_* cx_* i.year if collgrad==1 // 大学组
est store No
*-列示分组估计结果
esttab Yes No, nogap mtitle(Yes_Coll No_Coll) ///
star(* 0.1 ** 0.05 *** 0.01) s(r2 N)
*-似无相关估计
suest Yes No
*-组间差异检验
test [Yes_mean]cx_ttl_exp = [No_mean]cx_ttl_exp
test [Yes_mean]cx_hours = [No_mean]cx_hours
然后,我们数一下在这 10000 个随机数中,有多个是大于 (我们实际观察到的数值),命令为 count if d<-0.018。一共有 4963 个观察值大于 ,由此可得,经验 p 值 = 4963/10000 = 0.4963。这与采用 normal() 函数得到的结果 (0.4928) 非常接近。
. preserve
. clear
. set obs 10000
. set seed 1357
. gen d = rnormal() // d~N(0,1) 服从标准正态分布的随机数
. sum d, detail
. count if d<-0.018
. dis "Empirical p-value = " 4963/10000
. restore
C. 经验样本 (empirical sample)
上例中,我们假设 服从标准正态分布,从而可以通过 Monte Carlo 模拟的方式产生 10000 个观察值,这事实上是构造了一个经验样本。但多数情况下,我们并不知道 的分布特征,此时无法使用 Monte Carlo 模拟。然而,若假设抽样样本 (sample) 是从母体 (population) 中随机抽取的,则可以通过抽样样本中二次抽样得到经验样本 (empirical sample),这些经验样本也可以视为对母体的随机抽样。
若抽样过程中为无放回抽样 (sampling with no replacement),则相应的检验方法称为 “组合检验(permutation test)”;若抽样过程中为有放回抽样 (也称为可重复抽样,sampling with replacement),则对应的检验方法称为 “基于 Bootstrap 的检验”。
* 数据处理
. sysuse "nlsw88.dta", clear
. gen agesq = age*age
. drop if race==3
. gen black = 2.race
. global xx "ttl_exp married south hours tenure age agesq"
*-检验
. bdiff, group(black) model(reg wage $xx) reps(1000) detail
相关说明如下:
选项 group() 中填写用于区分组别的类别变量(若有多个组,可以预先删除不参与比较的组,类似于上面的 drop if race==3 命令);
选项 model() 用于设定回归模型,即上面提到的模型 (3a) 或 (3b);
选项 reps(#) 用于设定抽样次数,即上文提到的 K ,通常设定为 1000-5000 次即可;
附加 detail 选项,可以进一步列表呈现两组的实际估计系数 和 ;
上述过程大约用时 13 秒,结果如下:
- Permutaion (1000 times)- Test of Group (black 0 v.s 1) coeficients difference
Variables | b0-b1 Freq p-value
-------------+-------------------------------
ttl_exp | 0.007 490 0.490
married | -0.824 920 0.080
south | 1.411 5 0.005
hours | 0.010 344 0.344
tenure | -0.006 512 0.488
age | -1.579 751 0.249
agesq | 0.022 218 0.218
_cons | 28.051 267 0.267
---------------------------------------------
可以看到,ttl_exp 的经验 p 值为 0.49,表明白人和黑人组的 ttl_exp 系数不存在显著差异;married 变量的 p 值为 0.08,我们可以在 10% 水平上拒绝原假设。细心的读者会发现,该变量对应的 Freq = 920,为什么?(答案在上面 Step 5 处)。
*-数据预处理(可以忽略)
sysuse "nlsw88.dta", clear
gen agesq = age*age
*-分组虚拟变量
drop if race==3
gen black = 2.race
*-删除缺漏值
global xx "ttl_exp married south hours tenure age* i.industry"
qui reg wage $xx i.race
keep if e(sample)
*-生成行业虚拟变量
drop if industry==2 // 白人组中没有处于 Mining (industry=2) 的观察值
tab industry, gen(d) //手动生成行业虚拟变量
local dumind "d2 d3 d4 d5 d6 d7 d8 d9 d10 d11" //行业虚拟变量
global xx "c.ttl_exp married south c.hours c.tenure c.age c.agesq `dumind'"
*-permutation test
bdiff, group(black) model(reg wage $xx) reps(1000) detail
F. 面板数据
若原始数据为面板数据,通常会采用 xtreg, xtabond 等考虑个体效应的方法进行估计。抽样过程必须考虑面板数据的特征。在执行 bdiff 命令之前,只需设定 xtset id year,声明数据为面板数据格式,则抽样时便会以 id (公司或省份代码) 为单位,以保持 id 内部的时序特征。
*-Panel Data (sample by cluster(id))
*-数据预处理
. webuse "nlswork.dta", clear
. xtset id year //声明为面板数据,否则视为截面数据
. gen agesq = age*age
. drop if race==3
. gen black = 2.race
*-检验
. global x "ttl_exp hours tenure south age agesq"
. local m "xtreg ln_wage $x, fe" //模型设定
. bdiff, group(black) model(`m') reps(1000) bs first detail