论文解读与复现｜排污权交易机制是否提高了企业全要素生产率？

本文及相关材料由会员 @十月 投稿。

原文信息：任胜钢,郑晶晶,刘东华,陈晓红.排污权交易机制是否提高了企业全要素生产率——来自中国上市公司的证据[J].中国工业经济,2019(05):5-23.

导读

SO2排放是大气污染的主要来源，早期我国政府主要采用命令性手段来控制排放。1987年颁布了《大气污染防治法》，2002年环保部发布规定了电厂锅炉烟气脱硫的技术路线，但一刀切的环境规制手段限制了企业的生产积极性和治污的灵活性。

2000年我国开始采用市场化手段进行环境规制，2002年，国家环保部在山东、山西、江苏、河南4省，上海、天津、柳州3市，以及华能集团开展SO2排污有偿使用和交易政策的探索，即“4+3+1”项目。2007年正式在江苏、天津、浙江、河北、山西、重庆、湖北、陕西、内蒙古、湖南、河南11个省份开展试点工作。2017年SO2排污权交易试点在全国铺开。

波特假说认为环境规制会带来企业生产率的提高。排污权交易一方面可以促使企业通过技术创新提高全要素生产率，另一方面可以促使企业通过优化资源配置效率来提高全要素生产率。

文章以 2007 年 SO2 排放权交易试点政策为准自然实验，研究环境规制对企业全要素生产率的影响，以检验“波特假说”在中国是否成立。首先在省份和城市层面上运用双重差分(DID)法检验SO2排放权交易政策是否有效降低了试点地区的SO2排放。其次，在上市公司层面上运用DID模型检验SO2排放权交易政策对企业全要素生产率的影响，并进行了平行趋势检验。文章进一步在试点和非试点地区进行三重差分剔除其他排污权政策的影响，并运用更换全要素生产率测算方法、安慰剂检验、工具变量、剔除企业逃离行为的方法进行稳健性检验。最后，文章进行了机制分析和异质性分析。

样本和数据

样本为2004-2015年A股上市公司中的SO2排放企业。企业数据来自CSMAR数据库，企业专利数据来自佰腾网。

区分排放行业：以《污染防治规划》（环发[2008]1号）中规定的重点监测行业，电力、化工、石化、钢铁、有色、医药、造纸、冶金、水泥、玻璃等为SO2排放行业，其余为非SO2排放行业。区分排放企业：以企业年报和企业社会责任报告CSR中的SO2排放信息来确定。

DID样本为296家SO2排放企业，试点和非试点分别为111和185家。DDD样本为非SO2排放行业上市公司，共246家。

全要素生产率。TFP 测算采用Levinsohn and Petrin(2003)提出的半参数方法。估计代码及绘图：

use 基准回归数据.dta, clear
sum
xtset company year
gen lnzjz=ln(zjz+1)
gen lnlabor=ln(labor+1)
gen lnzjtr=ln(zjtr+1)
gen lncapital=ln(capital1+1)
gen tt=time*treat
gen lnzc=ln(zc+1)
gen lnzlb=ln(zlb+1)
gen lnaj=ln(aj+1)
gen lnrev=ln(revenue+1)

*******************************************
/* 图1:试点和非试点地区的平均TFP */
* 安装 levpet: net install st0060.pkg, from(http://www.stata-journal.com/software/sj4-2/)
levpet lnzjz, free(lnlabor) proxy(lnzjtr) capital(lncapital)
predict tfp,omega
gen lntfp=ln(tfp)
bys year treat: egen tfpmean = mean(lntfp)
form tfpmean %6.2f
gr tw ///
	conn tfpmean year if treat == 1, lp(solid) msym(o) mlabpos(6) || ///
	conn tfpmean year if treat == 0, lp(dash) msym(o) mlabpos(12) ///
	xla(2004(1)2015) yla(5.6(.2)6.6) xti("") ///
	leg(order(1 "试点地区" 2 "非试点地区") pos(10) ring(0)) ///
	xsize(20) ysize(12) scale(1.4)

图

实证结果

省级和城市层面的政策有效性检验

作者首先构建DID模型检验了SO2排放权交易政策对省份和地级市的SO2排放量的影响：

省级或城市处于试点地区则 ，否则为。年份大于等于年，否则为。的系数捕获试点地区在政策期间相对于非试点地区的排放平均变化。显著为负表明政策有效。为控制变量，和为固定效应。

回归代码及表2结果输出：

/* 表2 排污权交易政策对地区SO2排放和工业GDP的影响 */
/* DID：省份 */
use 省份政策有效性数据.dta,clear
gen tt=time*treat
reg so2 tt labour kx xf fdi gdzc i.year i.province, cluster(province)
est store m1

reg gygdp tt labour kx xf fdi gdzc i.year i.province, cluster(province)
est store m2

esttab m1 m2 using table2.rtf, b(%6.4f) se(%6.4f) r2(%9.4f) nogap compress star(* 0.1 ** 0.05 *** 0.01) indicate("RegionFE=*.province" "YearFE=*.year") scalar(N) title(table.2 col.1&3) replace

/* DID：城市 */
use 地级市政策有效性数据.dta,clear
gen tt = time*treat
reg so2 tt labour kx xf fdi gdzc i.year i.city, cluster(city)
est store m1

reg gygdp tt labour kx xf fdi gdzc i.year i.city, cluster(city)
est store m2

esttab m1 m2 using table2.rtf, b(%6.4f) se(%6.4f) r2(%9.4f) nogap compress star(* 0.1 ** 0.05 *** 0.01) indicate("RegionFE=*.city" "YearFE=*.year") scalar(N) title(table.2 col.2&4) append

表

排污权交易政策对企业全要素生产率的影响

同样，为了检验政策对企业全要素生产率的影响，构建如下DID模型：

为企业全要素生产率，表示企业位于试点地区，为行业固定效应，其他含义同前。捕捉了排污权交易政策对企业全要素生产率的影响，预期应显著为正。 回归代码及表 3 结果输出：

/* 表3 SO2试点对SO2排放企业的TFP */
set matsize 5000

use 基准回归数据.dta, clear
sum
xtset company year
gen lnzjz=ln(zjz+1)
gen lnlabor=ln(labor+1)
gen lnzjtr=ln(zjtr+1)
gen lncapital=ln(capital1+1)
gen tt=time*treat
gen lnzc=ln(zc+1)
gen lnzlb=ln(zlb+1)
gen lnaj=ln(aj+1)
gen lnrev=ln(revenue+1)

levpet lnzjz, free(lnlabor) proxy(lnzjtr) capital(lncapital)
predict tfp, omega
gen lntfp=ln(tfp)
bys year treat: egen tfpmean = mean(lntfp)
form tfpmean %6.2f

reg lntfp tt i.year i.area i.ind if so2==1,cluster(area)
est store m1

reg lntfp tt zcsy lf age owner sczy lnaj lnlabor lnzlb i.year i.area i.ind if so2==1,cluster(area)
est store m2

xtreg lntfp tt i.year i.company if so2==1, cluster(area)
est store m3

xtreg lntfp tt zcsy lf age owner sczy lnaj lnlabor lnzlb i.year i.company if so2==1, cluster(area)
est store m4

esttab m1 m2 m3 m4 using table3.rtf, b(%6.4f) se(%6.4f) r2(%6.4f) nogap compress star(* 0.1 ** 0.05 *** 0.01) indicate("Year FE =*.year" "Region FE =*.area" "Ind FE =*.ind" "Id FE =*.company") scalar(N) replace

表

平行趋势检验

使用双重差分的前提是处理组和对照组满足平行趋势假设，即在没有政策干预的情况下，结果变量在处理组和对照组的发展趋势一致。基准回归结果反映的是政策实施对企业全要素生产率的平均影响，没有反映政策在不同时间的影响差异。参考 Jacobson et al.(1993)的时间研究法对试点政策的动态效应进行实证：

对上述模型进行估计，并将 的系数和 95% 置信区间上下界绘制动态效应图。代码和图 2 绘制如下：

/* 平行趋势检验，2007年为基准年 */
gen tt1=treat*year04
gen tt2=treat*year05
gen tt3=treat*year06
gen tt5=treat*year08
gen tt6=treat*year09
gen tt7=treat*year10
gen tt8=treat*year11
gen tt9=treat*year12
gen tt10=treat*year13
gen tt11=treat*year14
gen tt12=treat*year15

reg lntfp tt1 tt2 tt3 tt5 tt6 tt7 tt8 tt9 tt10 tt11 tt12 zcsy lf age owner sczy lnaj lnlabor lnzlb i.year i.area i.ind if so2==1,cluster(area)
esttab, b(%6.4f) t(%6.4f) ar2(%9.4f) nogap compress star(* 0.1 ** 0.05 *** 0.01) drop(*.year *.area *.ind) scalar(N)

/* 95%CI */
matrix A=J(12,7,0)
forvalues i=1/12{
	if `i' != 4 {
	local k=2003+`i'
	scalar a`i'=_b[tt`i']
	scalar b`i'=(_b[tt`i']+invttail(e(df_r),0.025)*_se[tt`i'])
	scalar c`i'=(_b[tt`i']-invttail(e(df_r),0.025)*_se[tt`i'])
	mat A[`i',1]=`k'
	mat A[`i',2]=a`i'
	mat A[`i',3]=b`i'
	mat A[`i',4]=c`i'
	}
   else {
   	mat A[`i',1]=2007
   }
}

xtreg lntfp tt1 tt2 tt3 tt5 tt6 tt7 tt8 tt9 tt10 tt11 tt12 zcsy lf age owner sczy lnaj lnlabor lnzlb i.year if so2==1, fe cluster(area)

forvalues i=1/12{
	if `i' != 4 {
	scalar a`i'=_b[tt`i']
	scalar b`i'=(_b[tt`i']+invttail(e(df_r),0.025)*_se[tt`i'])
	scalar c`i'=(_b[tt`i']-invttail(e(df_r),0.025)*_se[tt`i'])
	mat A[`i',5]=a`i'
	mat A[`i',6]=b`i'
	mat A[`i',7]=c`i'
	}
}
mat l A

/* 图2 SO2排污权政策的动态效应 */
mat2txt, matrix(A) saving(DID_coe_ci) replace
insheet using DID_coe_ci.txt,clear
keep c1-c7
ren c1 year
* renvarlab 安装：ssc install renvarlab
renvarlab c2-c4 / coef up low
renvarlab coef up low, postfix(A)
renvarlab c5-c7 / coef up low
renvarlab coef up low, postfix(B)

gr two (sc coefA year, msymbol(O)) ///
	(rcap upA lowA year if year != 2007), ///
	xline(2007, lc(black)lp(dash)) yline(0, lc(black) ///
		lp(solid) lw(medthick)) ///
	xlabel(2004(1)2015) xtitle("")  ///
	yla(0 "0" -0.2(.2)0.6) ///
	leg(order(1 "estimate" 2 "95%CI") pos(6) cols(2)) ///
	name(a, replace) tit("(a) OLS")
gr two (sc coefB year, msymbol(O)) ///
	(rcap upB lowB year if year != 2007), ///
	xline(2007, lc(black)lp(dash)) yline(0, lc(black) ///
		lp(solid) lw(medthick)) ///
	xlabel(2004(1)2015) xtitle("") ///
	yla(0 "0" -0.2(.2)0.6) ///
	leg(order(1 "estimate" 2 "95%CI") pos(6) cols(2)) ///
	name(b, replace) tit("(b) FE")
gr combine a b, ysize(8) xsize(20) iscale(*1.6)

图

回归结果显示 tt1 到 tt6 的系数不显著，tt7 到 tt12 的系数显著，(a)(b)两个图中 2009 到 2010 年的系数发生了明显的跃升，说明在政策实施前处理组和对照组之间不存在显著差异，满足平行趋势假设，并且政策的影响存在 2 年的滞后效应。

三重差分

为了排除其他排污权政策可能带来的影响，作者选择制造业和采矿业中属于非 排放行业的企业为三重差分的另一对处理组和对照组。非SO2排放行业的企业不受到 排污权政策的影响，因此这一对处理组和对照组的差异只会来自于其他政策的影响。将之前的处理组和对照组的差异（SO2排放企业）减去这一对处理组和对照组的差异（非SO2排放企业），得到 排污权交易政策的净效应。三重差分模型如下：

其中， 为分类变量，样本企业属于 排放行业时 ，否则为 。 的系数 表示 排污权交易政策的净效应，系数为正说明促进了企业全要素生产率。回归代码如下：

/* 三重差分 */
use 基准回归数据.dta, clear
sum
xtset company year
gen lnzjz=ln(zjz+1)
gen lnlabor=ln(labor+1)
gen lnzjtr=ln(zjtr+1)
gen lncapital=ln(capital1+1)
gen tt=time*treat
gen lnzc=ln(zc+1)
gen lnzlb=ln(zlb+1)
gen lnaj=ln(aj+1)
gen lnrev=ln(revenue+1)

levpet lnzjz, free(lnlabor) proxy(lnzjtr) capital(lncapital)
predict tfp, omega
gen lntfp=ln(tfp)
bys year treat: egen tfpmean = mean(lntfp)
form tfpmean %6.2f

gen ttt=tt*so2
gen times=time*so2
gen treats=treat*so2

reg lntfp ttt tt treats times so2 i.year i.area i.ind, cluster(area)
est store m1
reg lntfp ttt tt treats times so2 zcsy lf owner age sczy lnaj lnlabor lnzlb i.year i.area i.ind ,cluster(area)
est store m2
xtreg lntfp ttt tt treats times so2 i.year i.company, cluster(area)
est store m3
xtreg lntfp ttt tt treats times so2 zcsy lf owner age sczy lnaj lnlabor lnzlb i.year i.company, cluster(area)
est store m4
/* 表4 三重差分结果 */
esttab m1 m2 m3 m4, b(%6.4f) se(%6.4f) r2(%9.4f) nogap compress star(* 0.1 ** 0.05 *** 0.01) indicate("时间效应 =*.year" "地区效应 =*.area" "行业效应 =*.ind" "个体效应 =*.company") scalar(N) replace

稳健性检验

GMM测算全要素生产率

Blundell and Bond(1998)通过GMM估计加入工具变量来解决模型中的内生性问题，文章采用GMM法估计企业全要素生产率，估计代码如下：

/* GMM测算全要素生产率 */
xtivreg2 lnzjz lnlabor (lncapital= l.lncapital) ,gmm2s fe robust
// 回归得到两个系数0.3324、0.5646
gen lntfp1=lnzjz-0.3324*lnlabor-0.5646*lncapital
// 得到GMM测算的tfp1
reg lntfp1 tt i.year i.area i.ind if so2==1,cluster(area)
est store m1

reg lntfp1 tt zcsy lf age owner sczy lnaj lnlabor lnzlb i.year i.area i.ind if so2==1,cluster(area)
est store m2

xtreg lntfp1 tt i.year i.company if so2==1, cluster(area)
est store m3

xtreg lntfp1 tt zcsy lf age owner sczy lnaj lnlabor lnzlb i.year i.company if so2==1, cluster(area)
est store m4
/* 表5 GMM 测算全要素生产率 */
esttab m1 m2 m3 m4, b(%6.4f) se(%6.4f) r2(%9.4f) nogap compress star(* 0.1 ** 0.05 *** 0.01) indicate("时间效应 =*.year" "地区效应 =*.area" "行业效应 =*.ind" "个体效应 =*.company") scalar(N)
/* 保存当前样本数据，方便安慰剂抽样时取用 */
save data.dta, replace

安慰剂检验

作者通过随机分配试点省份进行安慰剂测试，从31个省份中随机抽取11个省份作为处理组，其他作为对照组，重复500次抽样，并进行基准回归。表3第2列的真实估计（0.2708，对应图3中的红线）为安慰剂测试中的异常值时，说明文章估计结果不是由省份-行业-年份中的不可观测因素驱动的。抽样代码及安慰剂检验绘图如下：

/* 安慰剂检验 */
use data, clear
mat b = J(500,1,0)
mat se = J(500,1,0)
mat p = J(500,1,0)

* 500次随机抽样
forvalues i=1/500{
	di "`i'"
	qui {
		use data.dta, clear
		keep area
		duplicates drop area,force
		sample 11, count // 31个省份中随机抽取11个省份
		mkmat area, matrix(samplearea) //向量转化为矩阵
		use data.dta,clear
		xtset company year
		replace treat = 0
	  * 生成tt处理变量
		foreach j of numlist 1/11 {
			replace treat = 1 if (area == samplearea[`j',1])
			replace tt = time * treat
		}
		qui xtreg lntfp tt zcsy lf age owner sczy lnaj lnlabor lnzlb i.year if so2==1,fe cluster(area)
	  * 存储并计算所需回归结果
		mat b[`i',1] = _b[tt]
		mat se[`i',1] = _se[tt]
		mat p[`i',1] = 2*ttail(e(df_r),abs(_b[tt]/_se[tt]))
	}
}

/* 图3 安慰剂检验 */
svmat b, names(estimate)
svmat se, names(se)
svmat p, names(pvalue)

drop if pvalue1 == .

twoway ///
	kdensity estimate1|| ///
	scatter pvalue1 estimate1, ///
	msymbol(smcircle_hollow) ytitle(p-value) xtitle(estimates) ///
	legend(order(1 "kdensity of estimates" 2 "p-value") ring(0) position(10) cols(1)) ///
	yline(0.01) xline(0.2708,lc(red)lw(medium)) ///
	xlabel(-.4(.1).35) scale(1.2)

安慰剂检验

工具变量

依据 Hering and Poncet(2004)，作者使用通风系数作为工具变量，通风系数为风速和混合层高度的乘积，值越大意味着污染物扩散越快。通风系数越大，被选为试点的概率越小，因此满足相关性假设。同时通风系数由气象因素和地理因素决定，满足外生性假设。工具变量回归代码如下：

use data.dta, clear

/* 工具变量 */
gen tf1=tf*time
gen tf2=tf1*tf1
// 第二阶段的结果
ivregress 2sls lntfp (tt=tf1) zcsy lf age owner sczy lnaj lnlabor lnzlb i.year i.area i.ind if so2==1,cluster(area)
est sto m1

xtivreg lntfp (tt=tf1) zcsy lf age owner sczy lnaj lnlabor lnzlb i.year i.company if so2==1, vce(cluster area)
est sto m2
/* 表6 工具变量回归结果*/
esttab m1 m2, b(%6.4f) se(%6.4f) r2(%9.4f) ar2(%9.4f) nogap compress star(* 0.1 ** 0.05 *** 0.01) drop(*.year *.company *.ind *.area) scalar(N)

企业逃离

为了规避环境规制，一些企业可能从试点搬离到非试点地区。因此，作者删除了注册地址发生变化的企业，重新回归，结果依然稳健。代码和结果如下：

/* 企业逃离 */
reg lntfp tt i.year i.area i.ind if so2==1 & dzbq==0, cluster(area)
est sto m1

reg lntfp tt zcsy lf age owner sczy lnaj lnlabor lnzlb i.year i.area i.ind if so2==1&dzbq==0,cluster(area)
est sto m2

xtreg lntfp tt i.year i.company if so2==1&dzbq==0,cluster(area)
est sto m3

xtreg lntfp tt zcsy lf age owner sczy lnaj lnlabor lnzlb i.year i.company if so2==1&dzbq==0,cluster(area)
est sto m4
/* 表7 剔除搬迁的企业 */
esttab m1 m2 m3 m4, b(%6.4f) se(%6.4f) r2(%9.4f) nogap compress star(* 0.1 ** 0.05 *** 0.01) indicate("时间效应 =*.year" "地区效应 =*.area" "行业效应 =*.ind" "个体效应 =*.company") scalar(N)

*> --------------------------------------------------------------
*>                  (1)          (2)          (3)          (4)
*>                lntfp        lntfp        lntfp        lntfp
*> --------------------------------------------------------------
*> tt            0.2508***    0.2652***    0.2508***    0.2577***
*>             (0.0827)     (0.0707)     (0.0856)     (0.0755)
*> zcsy                       0.0132***                 0.0093***
*>                          (0.0022)                  (0.0017)
*> lf                        -0.0112*                  -0.0103*
*>                          (0.0057)                  (0.0057)
*> age                       -0.0025                    0.0043***
*>                          (0.0031)                  (0.0016)
*> owner                      0.0799                    0.3744*
*>                          (0.0649)                  (0.1926)
*> sczy                       0.0215***                 0.0059**
*>                          (0.0074)                  (0.0025)
*> lnaj                       0.0380                    0.0605**
*>                          (0.0333)                  (0.0306)
*> lnlabor                    0.1628***                 0.0009
*>                          (0.0360)                  (0.0563)
*> lnzlb                      0.0592***                -0.0150
*>                          (0.0167)                  (0.0221)
*> _cons         6.6344***    4.5657***    5.9667***    5.5387***
*>             (0.2407)     (0.4952)     (0.0501)     (0.4242)
*> 时间效应       Yes          Yes          Yes          Yes
*> 地区效应       Yes          Yes           No           No
*> 行业效应       Yes          Yes           No           No
*> 个体效应        No           No          Yes          Yes
*> --------------------------------------------------------------
*> N               3312         3311         3312         3311
*> R-sq          0.2330       0.3794
*> --------------------------------------------------------------
*> Standard errors in parentheses
*> * p<0.1, ** p<0.05, *** p<0.01

机制分析

全要素生产率是一个综合的索洛余值，可以分解为技术进步效应和资源配置效应两个方面。因此，文章从技术创新和资源配置两方面进行机制分析。

技术创新

作者选取专利申请数量作为技术创新的代理变量，检验试点政策对技术创新的影响。作者披露的数据中缺少 xfs 变量，此处复现的是剔除了 xfs 变量的结果，不影响文章结论，代码如下：

/* 传导机制 */

/* 技术创新 */
reg zlb tt zcsy zc lf iso lnlabor  sch i.year i.area i.ind if so2==1,cluster(area)
est sto m1
reg zlb tt zcsy zc lf iso lnlabor sch i.year i.company if so2==1,cluster(area)
est sto m2
/* 表8 列1 列2 */
esttab m1 m2, b(%6.4f) se(%6.4f) r2(%9.4f) nogap compress star(* 0.1 ** 0.05 *** 0.01) indicate("时间效应 =*.year" "地区效应 =*.area" "行业效应 =*.ind" "个体效应 =*.company") scalar(N)

*> ------------------------------------
*>                  (1)          (2)
*>                  zlb          zlb
*> ------------------------------------
*> tt           10.0728***   10.1961***
*>             (3.0683)     (3.0969)
*> zcsy          0.0059      -0.0226
*>             (0.0240)     (0.0161)
*> zc            0.0000**     0.0000*
*>             (0.0000)     (0.0000)
*> lf            0.3391       0.3386
*>             (0.2031)     (0.2110)
*> iso          11.7700       9.7404*
*>             (7.1682)     (5.5862)
*> lnlabor       0.5201       0.0945
*>             (0.9245)     (1.6389)
*> sch          -1.9631      -1.9186
*>             (2.4058)     (2.4490)
*> _cons         9.3578      25.0056
*>            (17.9859)    (30.1175)
*> 时间效应       Yes          Yes
*> 地区效应       Yes           No
*> 行业效应       Yes           No
*> 个体效应        No          Yes
*> ------------------------------------
*> N               3480         3480
*> R-sq          0.2737       0.5777
*> ------------------------------------
*> Standard errors in parentheses
*> * p<0.1, ** p<0.05, *** p<0.01

资源配置效率

作者采用两种方式验证资源配置效率，（1）借鉴徐升艳等（2018），在模型中控制技术创新，将全要素生产率中的技术进步效应剔除。（2）借鉴方军雄（2007）投资-投资机会”敏感性模型，利用资本配置效率替代资源配置效率。

/* 资源配置效率 */
gen ttroa=tt*roa
// 控制技术创新
reg lntfp tt zcsy lf age owner sczy lnaj lnlabor lnzlb i.year i.area i.ind if so2==1,cluster(area)
est sto m1

xtreg lntfp tt zcsy lf age owner sczy lnaj lnlabor lnzlb i.year i.area i.ind i.company if so2==1,cluster(area)
est sto m2
// 使用资本配置效率替代
reg tz ttroa roa tt age zlb lf sczy iso i.year i.area i.ind if so2==1,cluster(area)
est sto m3

xtreg tz ttroa roa tt age zlb lf sczy iso i.year i.company if so2==1,cluster(area)
est sto m4
/* 表8 列3 列4 列5 列6 */
esttab m1 m2 m3 m4, b(%6.4f) se(%6.4f) r2(%9.4f) nogap compress star(* 0.1 ** 0.05 *** 0.01) indicate("时间效应 =*.year" "地区效应 =*.area" "行业效应 =*.ind" "个体效应 =*.company") scalar(N)

异质性分析

文章从企业外部特征和内部特征两方面进行异质性讨论。具体而言，外部特征为环境执法力度，内部特征为企业所有制。作者分别汇报了国有和非国有样本的估计值，国有企业的估计系数小于非国有企业。环境执法强度高的估计系数大于环境执法强度低的系数。估计代码和结果如下：

/* 异质性检验 */
egen ajj = cut(aj), group(2)

// 国有和非国有
reg lntfp tt zcsy lf age sczy lnaj lnlabor lnzlb i.year i.company if so2==1&owner==0,cluster(area)
est sto m1

reg lntfp tt zcsy lf age sczy lnaj lnlabor lnzlb i.year i.company if so2==1&owner==1,cluster(area)
est sto m2

// 环境执法强度
reg lntfp tt zcsy lf age sczy owner lnlabor lnzlb i.year i.company if so2==1&ajj==0,cluster(area)
est sto m3

reg lntfp tt zcsy lf age sczy owner lnlabor lnzlb i.year i.company if so2==1&ajj==1,cluster(area)
est sto m4

/* 表9 异质性分析 */
esttab m1 m2 m3 m4, b(%6.4f) se(%6.4f) r2(%9.4f) nogap compress star(* 0.1 ** 0.05 *** 0.01) indicate("时间效应 =*.year" "个体效应 =*.company") scalar(N)

*> --------------------------------------------------------------
*>                  (1)          (2)          (3)          (4)
*>                lntfp        lntfp        lntfp        lntfp
*> --------------------------------------------------------------
*> tt            0.2950*      0.2471**     0.3720***    0.1282***
*>             (0.1539)     (0.1008)     (0.1285)     (0.0418)
*> zcsy          0.0066**     0.0112***    0.0101***    0.0089***
*>             (0.0025)     (0.0025)     (0.0020)     (0.0022)
*> lf            0.0054      -0.0101      -0.0048      -0.0027
*>             (0.0088)     (0.0068)     (0.0102)     (0.0061)
*> age          -0.0021       0.0979***   -0.0198       0.0237***
*>             (0.0157)     (0.0279)     (0.0120)     (0.0023)
*> sczy         -0.0088       0.0054*      0.0066***    0.0101
*>             (0.0204)     (0.0031)     (0.0023)     (0.0161)
*> lnaj          0.0789*      0.0588
*>             (0.0421)     (0.0353)
*> lnlabor       0.0963      -0.0371       0.0679      -0.0388
*>             (0.0855)     (0.0694)     (0.0725)     (0.0609)
*> lnzlb         0.0127      -0.0269      -0.0138      -0.0327
*>             (0.0338)     (0.0297)     (0.0350)     (0.0238)
*> owner                                  -0.7846***    0.6776***
*>                                       (0.1419)     (0.0380)
*> _cons         4.6504***    5.9384***    5.5502***    6.1183***
*>             (0.6430)     (0.7902)     (0.4545)     (0.5288)
*> 时间效应       Yes          Yes          Yes          Yes
*> 个体效应       Yes          Yes          Yes          Yes
*> --------------------------------------------------------------
*> N               1379         2100         1824         1655
*> R-sq          0.5496       0.5883       0.5667       0.6372
*> --------------------------------------------------------------
*> Standard errors in parentheses
*> * p<0.1, ** p<0.05, *** p<0.01

点击这里跳转到 RStata 短书平台获取附件：论文解读与复现｜排污权交易机制是否提高了企业全要素生产率？