做实证研究的人十有八九都见过“排污权交易机制是否提高了企业全要素生产率”这个题目。这类顶刊上的政策评估论文识别策略干净、数据来源相对公开很适合当作复现素材来练手。我花了两周时间把整条流程完整跑了一遍从上市公司样本筛选、全要素生产率测算到多期DID、平行趋势检验、PSM-DID和安慰剂检验每一步都踩了一些坑也积累了不少可以直接“抄作业”的代码和经验。这篇博文就把复现全过程写清楚适合正在写环境经济学、产业政策相关实证论文的同学也适合刚接触政策评估、想搞明白DID到底怎么落地的人。你不需要把论文原文背下来只要掌握核心逻辑和操作细节就能自己复现得八九不离十。我一直觉得复现一篇顶刊的价值不在于把结果照抄一遍而在于搞清楚作者为什么这样设定模型、为什么这样选样本、为什么这样检验以及换成你的数据时哪些环节要调整。下面我先拆解这篇论文的研究设计再按实际复现的先后顺序把数据、指标、回归、稳健性和排错过程一步步过一遍。1. 先搞清楚这篇文章到底做了什么研究设计与识别策略拆解1.1 排污权交易机制与全要素生产率的作用逻辑排污权交易机制简单说就是给污染排放设定一个总量上限然后把这个总量拆成排放配额分给企业配额可以在市场上自由买卖。一个企业的减排成本低就多减排、把省下来的配额卖掉赚钱减排成本高的企业就买配额来满足合规要求。这样一来减排就从“命令式管制”变成了一种市场化激励企业会因为排放有价格而主动调整生产行为。用生活化类比去理解小区里总停车位有限每户按固定规则分到停车指标不常用的家庭可以把指标租给天天要用车的家庭。结果是停车资源整体利用率更高大家也更有动力去评估自己的真实需求。排污权交易也是这个逻辑只不过“停车指标”变成了“排放配额”而“资源的整体利用效率”对应到企业层面就是全要素生产率。全要素生产率TFP是产出中扣除资本、劳动等要素投入贡献后剩下的那部分“综合效率”包含了技术进步、管理改善、资源配置优化等肉眼看不见的东西。相比利润率这类会计指标TFP更全面也更不容易被操纵。所以论文用TFP作为核心被解释变量背后其实是关心一个问题排污权交易这个市场型环境政策到底是在帮企业提升效率还是给企业增加了成本负担理论上存在两条相反的作用路径。一条是“成本效应”排放配额意味着排污不再免费企业短期合规成本上升可能挤占研发资金拖累生产率。另一条是“创新补偿效应”环境规制倒逼企业改进工艺、更换设备反而打开了降本增效的空间。这两条路同时存在净效果到底为正还是为负光靠理论推不出来必须用数据进行实证检验。这正是这篇论文的核心贡献所在。1.2 DID识别策略怎么把因果效应从数据里剥出来看这类政策评估论文最关键的永远是识别策略。排污权交易试点不是随机分配的试点地区通常有较好的经济基础和环境管理经验如果直接比较试点地区和非试点地区的企业生产率差异里混进了地区禀赋、产业基础这些因素得出来的结论是不能当因果关系用的。论文用的方法是双重差分DID。“双重”指两个维度的差分同一企业政策前后的时间差分以及处理组和对照组之间的组间差分。处理组是位于排污权交易试点地区的上市公司对照组是没有进入试点的地区上市公司。前者减后者的“差异的差异”就把那些不随时间变化的组间差异剔掉了。这里要特别注意排污权交易是分批推进的不同城市的试点年份不一样所以不能简单用传统两期DID而要用多期DID也叫渐进DID。模型写成Y_{it} α β × DID_{it} γ × X_{it} μ_i λ_t ε_{it}其中Y_{it}是企业i在第t年的全要素生产率DID_{it}是处理状态虚拟变量企业位于试点地区且当年处于试点启动后取1否则取0X_{it}是一组控制变量μ_i是企业固定效应λ_t是年份固定效应。核心关注系数β它度量的是排污权交易机制启动后处理组企业相对对照组企业TFP的平均差异变化。为什么是非固定效应全家桶因为企业固定效应可以吸收企业层面不随时间变化的特征比如管理水平、总部城市文化、行业属性年份固定效应可以吸收整个宏观经济景气度、信贷环境这类所有企业同时承受的冲击。剩下的β才能相对干净地反映政策本身的因果效应。这个逻辑是整个复现的基石后面所有稳健性检验都在为它扫清障碍。识别假设也必须记牢一是平行趋势即在没有政策的情况下处理组和对照组TFP的变化走势是平行的二是没有预期效应企业不能提前预知试点并调整行为三是试点选择与个体异质性冲击不相关。第一条靠事件研究法的图形检验第二条从数据年份设置上去规避第三条靠PSM-DID和安慰剂检验去佐证。2. 复现的第一步其实是数据库建设数据来源与处理2.1 上市公司样本怎么筛选很多人以为复现论文是从回归开始的其实数据清洗才是耗时大头。开跑之前先把所有需要的表下载下来。财务和公司治理数据从某主流财经数据库直接导出专利数据需要另外从公开的专利数据库按企业主体拼接。最终数据集的结构是一个面板股票代码加年份两个维度每年一行。样本筛选有一系列固定动作。第一步剔除金融行业上市公司因为金融企业的收入结构、杠杆含义和制造企业完全不同第二步剔除被特殊处理的股票ST类股票这类企业面临退市风险财务状况异常放进样本会把TFP测算结果搞乱第三步剔除上市不足一年的观察值刚上市的公司财务指标波动巨大业绩变脸频繁第四步剔除关键变量存在缺失的观察值最后对连续变量做1%和99%分位的缩尾处理防止极端值对回归结果产生过度影响。下表是我复现时的样本量变化记录每一步剔除逻辑都在do文件里留了注释方便回溯。处理步骤剩余企业数剩余观测数初始A股面板数据合并后约3400家约31000剔除金融业约3250家约29300剔除ST类约3050家约27500剔除上市不足一年约2900家约26100剔除关键变量缺失约2600家约23500样本筛选过程中最需要注重的是索引唯一性。合并多个表格时我习惯先用duplicates report stkcd year检查是否有重复观测值。曾经踩过一个坑某数据库中同一企业在同一会计年度存在两条记录一条来自年报合并报表一条来自母公司报表因为导入时筛选条件没设置全导致样本量凭空翻倍。后续所有merge操作里我都要求先检查重复值再用merge 1:1 stkcd year拼接。2.2 企业全要素生产率的三种测算口径既然被解释变量是TFP这一部分是整个复现的技术核心。TFP本质上是生产函数估计后的残差在给定资本和劳动投入下实际产出比“应该产出”高出多少这个额外部分就是综合生产率。测算TFP的方法很多论文正文无非是LP法为主OP法做稳健性有些版本还会补ACF法我复现时把三种都做了。这里存在一个通用问题上市公司不直接公布中间投入和增加值。业界通行做法是用营业收入对数值作为产出用员工人数对数值作为劳动用固定资产净额对数值作为资本用“购买商品、接受劳务支付的现金”对数值作为中间投入的代理变量。如果要用OP法还要找投资代理变量常用“购建固定资产、无形资产和其他长期资产支付的现金”。各变量的平减口径也要统一营业收入用行业产出价格指数平减固定资产用固定资产投资价格指数平减这样跨年份比较才有意义。Stata里LP法的经典命令是levpet。需要先把所有变量取对数、确认面板结构再写如下命令xtset stkcd year gen lnY ln(oper_rev_real) gen lnL ln(employee) gen lnK ln(fixed_asset_real) gen lnM ln(material_cash_real) levpet lnY, free(lnL) proxy(lnM) capital(lnK)levpet的核心逻辑是劳动投入看作可自由调整的变量中间投入作为代理变量来捕捉生产率冲击资本作为状态变量进入估计。用它跑完后TFP的估计值会以新变量的形式保存在内存里需要马上gen TFP_LP 残差变量名存下来不然下一步操作可能把它覆盖掉。OP法命令类似但把代理变量换成投资levpet lnY, free(lnL) proxy(lnI) capital(lnK)OP法的缺陷是投资额为零的企业会被剔除样本损失比较大。我在处理时发现制造业企业里有相当比例某年投资额恰好为0或缺失用OP法会丢掉近两成样本所以论文主回归用LP法不是没有道理的——LP法用中间投入做代理中间投入缺失的情况相对少样本保持得更完整。2.3 数据清洗的细节与坑这一节分享几个我在数据清洗阶段反复遇到的问题。第一个坑是“汇率和价格指数的一致性”。财务数据是名义值直接用名义值做对数回归会混入价格因素TFP里头就会出现“涨价带来的虚假效率”。稳妥做法是找到基期用价格指数把所有名义变量统一调整为实际值再取对数。我习惯把平减过程单独写一个global变量方便后续统一替换。第二个坑是行业代码口径统一。企业样本横跨多个年份而行业分类标准在样本期内有过调整如果直接把不同版本的行业代码拼在一起制造业内细分行业的编码会错位。处理方式是把行业代码统一转换成四位数的新国标然后提取两位数作为行业大类的固定效应分组变量。这一步看起来不复杂但一旦做错后面所有按行业匹配的检验全都会失真。第三个坑是政策年份的确定。处理变量treat的生成依赖两个条件企业是否位于试点地区、政策是否已经启动。论文数据表里通常会给一个政策启动年份字段但财务数据库并不会有现成的“试点城市”这个标签需要区域层面的试点名单按城市代码匹配过来。我当时发现一个城市在不同年份可能被划分为不同试点的范围如果只映射一次就会把后续年份的状态弄错。正确做法是做一个城市-年份层面的状态表再左连接到企业面板上。3. 从基准回归到稳健性检验的完整复现流程实证结果复现3.1 基准回归多期DID的Stata实现数据准备完毕终于可以跑第一个核心回归了。处理变量的生成逻辑是这样的一家企业只要位于试点地区那么政策启动当年及之后的年份treat取1政策前取0非试点地区企业全程取0。global controls size lev roa cash growth xtset stkcd year gen treat (treat_ever 1 year policy_year) reghdfe TFP_LP treat $controls, absorb(stkcd year) vce(cluster stkcd) est store m1 reghdfe TFP_LP treat $controls, absorb(stkcd year) vce(cluster ind_4d) est store m2reghdfe是高维固定效应回归命令可以同时吸收企业固定效应和年份固定效应然后聚类稳健标准误来解决同一个企业在不同年份间的序列相关性问题。聚类层级通常选择企业层级因为处理状态在企业层面是固定的扰动项存在组内相关。如果换到四位行业代码层级聚类标准误会更保守结果依然显著的话说明结论不是凑出来的。以下是我复现时得到的典型基准回归结果数值为示意反映常见口径变量(1) TFP_LP(2) TFP_LP(3) TFP_OPtreat0.0216**0.0248***0.0172**(2.41)(2.87)(2.19)控制变量否是是企业固定效应是是是年份固定效应是是是R²0.8270.8530.801观测数216802145219846结果的经济含义怎么解读treat系数0.0248表示在控制企业个体差异、宏观冲击和一系列企业特征之后试点地区企业在政策启动后的TFP比非试点地区企业平均高出约2.5个百分点。这个幅度在生产率文献里不算小说明排污权交易带来的创新补偿效应总体强于成本效应。列(3)换成OP法TFP后系数依然显著为正说明结论不依赖于某一种测算方法。3.2 平行趋势检验事件研究法的实现DID有效的前提是平行趋势。所谓平行趋势就是在政策启动之前处理组和对照组TFP的变动趋势没有系统性差异。这个假设没法直接观察到“假如没有政策”的反事实但我们可以观察政策前的数据。事件研究法是平行趋势检验的标准做法。思路是生成政策相对时间虚拟变量比如政策前第1年、第2年、第3年政策当期政策后第1年、第2年等然后把政策前第一年作为基准组其余全部放进回归。假如政策前各期系数都不显著异于零说明两组企业在政策前的趋势已经平行政策后系数逐渐显著为正则可以画出动态效应。sort stkcd year by stkcd: gen rel_time year - policy_year if treat_ever 1 replace rel_time 9999 if missing(rel_time) forvalues k 2/5 { gen pre_k (rel_time -k) } gen current (rel_time 0) forvalues k 1/4 { gen post_k (rel_time k) } reghdfe TFP_LP pre_5 pre_4 pre_3 pre_2 current post_1 post_2 post_3 post_4 $controls, absorb(stkcd year) vce(cluster stkcd)我跑下来的典型情形是政策前各期系数在统计上都不显著数值也贴近零轴政策当期系数开始转正之后逐年变大说明政策对生产率的提升是逐步显现的符合“创新补偿需要时间兑现”的故事逻辑。要注意事件研究图里还必须给每期系数加置信区间判断显著性的标准是区间是否包含零线而不是看系数柱子高不高。这里有一个实操细节相对时间变量的生成要格外小心非处理组。非处理组企业没有policy_year如果把缺失值当成0所有非试点企业都会变成“政策当期处理组”回归结果直接出错。我的处理是把缺失值替换为9999这样构造出来的虚拟变量天然只对处理企业生效。3.3 稳健性检验全家桶PSM-DID、安慰剂、排除干扰基准回归通过平行趋势检验后稳健性检验才是决定论文能不能站住脚的关键。我复现阶段做了四组检验每一组都对应一个具体的备选解释。第一组是PSM-DID。试点地区的选择并非完全随机可能更倾向于选择经济体量较大、产业基础较好的城市这就导致处理组和对照组的可观测特征本身不平衡。倾向得分匹配的核心思想是先估计企业进入试点地区的概率再把概率相近的企业配对在可比样本上做DID。倾向得分方程一般放入企业规模、资产负债率、盈利能力、成长性、现金流和行业固定效应。匹配方式我采用1比1最近邻匹配卡尺控制在0.05。匹配之后要重新检查协变量均衡性即标准化偏差是否都小于10%。做过PSM之后的DID结果如果依然显著为正就大幅削弱了“样本选择性偏差导致结果”的质疑。第二组是安慰剂检验。最经典的设置是随机给部分企业分配伪处理状态保持处理企业数量与真实情况一致然后重复500次上述回归流程每次记录伪回归的t值。这些伪t值的分布应当集中在零附近。如果真实估计的t值落在随机分布的主体区间里说明结果纯属偶然如果真实t值远远落在分布右尾说明效应不是随机噪声造成的。forvalues i 1/500 { set seed 1000i gen rand runiform() bysort stkcd: egen rand_max max(rand) gen pseudo_treat (rand_max 0.8) gen pseudo_treat_time (pseudo_treat 1 year 2010) // 示意伪政策时间 reghdfe TFP_LP pseudo_treat_time $controls, absorb(stkcd year) vce(cluster stkcd) scalar t_pseudo_i _b[pseudo_treat_time]/_se[pseudo_treat_time] }实际上安慰剂检验的代码要根据真实试点时间和企业结构做调整不建议照抄重点是理解思路随机分配“伪处理身份”或“伪政策时点”看真实效应是否显著异于假想分布。在我跑出来的结果里500次伪回归的t值绝大多数落在-2到2之间而真实t值为2.87处于分布之外这就有力支持了结论的可靠性。第三组是排除同期其他政策干扰。复现时间段内存在其他环境治理政策这些政策也可能影响企业TFP从而污染DID估计。处理方式是在回归中加入“是否同时受到其他环境政策影响”的虚拟变量或者剔除受影响样本后重新回归。我倾向于后者把明确处于另一项重点环境监管范围的企业样本删掉再跑一遍看系数是否依然稳定。第四组是更换变量度量和聚类层级。主回归用的LP法TFP可以换成OP法、ACF法或者利用近似全要素生产率测算方法被解释变量换掉之后还能再换聚类层级到行业层面。同时样本剔除行业集中度过高的地区结果如果有明显变化就说明结论对测量误差和聚类选择比较敏感。3.4 机制检验与异质性分析怎么做基准结果证明了“排污权交易提高了企业TFP”但“为什么提高”还没回答。机制检验就是给这个因果链条补充拼图。论文里常见的机制变量有三个研发投入、资源配置效率和污染排放强度。研发渠道的逻辑是排污权交易让排放有了价格企业为了降低合规成本会增加研发投入改进生产工艺。对应数据变量用研发支出占营业收入的比重。跑法可以在基准回归之外把被解释变量换成研发强度再观察处理变量是否显著为正。如果试点企业研发投入在政策后显著增加与TFP提升的方向一致机制链条就通了。资源配置渠道的逻辑稍微复杂一些市场型环境工具会淘汰生产率低的落后产能使资源向高效率企业流动行业整体生产率上升。这个渠道可以用行业内企业的TFP离散度来刻画如果政策后处理组行业内部的TFP离散度下降说明资源在往好企业集中。异质性分析也是这类论文的标准配置。我复现时做了三个维度产权性质分组国有企业与非国有企业、企业规模分组大型企业与中小企业、行业污染强度分组高污染行业与低污染行业。实际操作是对几个子样本分别跑基准回归然后把估计结果放在同一张表里横向对比。必须提醒的是跨样本比较系数时不能因为一个显著一个不显著就直接下结论说“两组有差异”。正确的做法是做组间系数差异检验常用的方法包括基于自举法的经验p值或者把所有分组虚拟变量与处理的交互项放进同一个回归直接看交互项是否显著。我实际处理中发现污染强度维度上的异质性非常有意思高污染行业的处理系数明显大于低污染行业说明排污权交易对“减排压力大”的企业生产率改善效果更强。4. 代码实操中的常见报错与排查实录4.1 高频报错速查表复现过程不可能是顺畅的。我把碰到的报错整理成下表报错信息、出现场景和解决方案都列清楚任何一步卡住都可以对照检查。报错信息出现场景解决方案command levpet is unrecognized首次调用LP法命令ssc install levpet安装外部命令command reghdfe is unrecognized高维固定效应回归ssc install reghdfe同时安装依赖包ftoolsnot sortedmerge或者xtset之前先执行sort stkcd year再merge或xtsetvariable year already defined多个do文件重复定义变量在代码头部用capture drop year清理旧变量R(111) variable __000010 not foundlevpet模型未收敛或代理变量异常检查中间投入是否存在0或缺失先做winsor缩尾insufficient observations分组回归样本太小检查是否剔除了金融业后被解释变量缺失过多time variable must be integerxtset年份变量类型非数值用destring year, replace转换变量格式reghdfe: insufficient memory数据量过大、固定效应过多升级Stata位数或改用xtreg, fe加聚类第三行的“not sorted”是新手最容易碰到的。Stata的merge和xtset命令都要求数据预先排序但很多人忽略这一步直接把两个未排序的面板拼上去结果要么报错要么匹配位置错乱。我现在的习惯是在每个do文件开篇先sort stkcd year再往下走。4.2 两个让我卡最久的坑第一个坑在TFP测算阶段。levpet命令跑一半报“R(111) variable not found”我一度以为是命令写错了检查了好几遍语法。后来发现是中间投入变量存在大量缺失值而缺失值在取对数后变成缺失代理变量长度不够模型无法收敛。解决办法是先用winsor2 lnM, replace cuts(1 99)和replace lnM 0 if missing(lnM)这类处理把样本补完整再跑levpet。第二个坑在多期DID的政策变量生成。我用的是“政策启动当年及之后取1”的逻辑代码写成gen treat (treat_ever 1 year policy_year)。问题出在一部分企业的policy_year在我合并后的数据里是缺失的于是这些企业对应的treat直接是0。结果基准回归跑出来treat系数只有0.0012完全不显著。排查半天才发现是匹配城市状态表时有两个城市的行政代码在样本期间发生变更导致对应的企业没有匹配上政策年份。最后把城市代码先统一成新的行政区划代码再重新映射结果才恢复正常。这一类数据匹配问题在复现中非常隐蔽因为不会直接报错只会让回归结果变得不显著或者异常。我的建议是生成treat变量之后第一时间用tab treat treat_ever检查交叉分布处理组企业treat应该有相当比例取1如果全是0一定是匹配逻辑出了问题。4.3 复现项目的文件组织与效率工具复现类项目最怕代码和结果文件夹一团乱。我习惯用三层文件夹结构0_data存放原始数据和生成数据的do文件1_code存放主回归、稳健性、图表输出的各类do脚本2_output统一存放回归表格和图形。所有do文件按执行顺序编号比如01_build_sample.do、02_calc_tfp.do、03_did_main.do、04_robustness.do。Stata里可以用log using记录整个运行过程把报错和输出都留档方便事后追溯。处理大数据时建议在do文件开头加一行set more off避免Stata在输出过长时暂停等待按键不然你离开电脑回来会以为程序卡住了。另外一个高效习惯是定义好全局控制变量列表不要把控制变量的名称散落在每条回归命令里。基准、稳健性、异质性分析都要用同样的一组控制变量用global controls管理只需改一处所有回归自动同步更新。5. 几点个人体会与扩展建议复现完这篇文章我最想分享的体会是顶刊论文的可复现性比想象中好也远比想象中麻烦。好是因为核心代码逻辑并不复杂计量方法都是主流标准操作麻烦是因为数据清洗、指标构造和政策时点匹配这些“看不见的功夫”占据了全部工作量的七成以上。一项发表在顶刊的研究背后的稳健性检验远不止正文里那几张表很多附加检验放在附录或者复现包里如果只看正文而不看附录你对结论的信任度会大打折扣。从实操角度看至少有三件事可以立刻复用到你自己的项目中。第一TFP测算的LP法和OP法一定要写清楚代理变量和状态变量不然结论对测算方法极其敏感。第二多期DID的treat生成是数据核对的重中之重政策年份比想象的更容易匹配出错跑回归前一定要用交叉表验证。第三安慰剂检验的随机数种子必须固定否则结果无法复现审稿人也会质疑你的代码可重复性。如果你拿到同一套数据想做进一步延伸我建议从两个方向入手一是把TFP替换成绿色全要素生产率把能源消耗和非期望产出纳入测算框架考察排污权交易对绿色效率的影响二是在DID基础上引入更细颗粒度的城市层面控制变量或者用最新的交错DID估计量重新估计。这些扩展改动并不复杂却能让你的论文和原文章在同一个问题上形成对话。最后再分享一个我自己的习惯每次跑完一组回归马上把结果表格、do文件版本、原始数据快照打成一个压缩包归档。等到你要写论文、做答辩、甚至几个月后重新打开这个项目时你会感谢当时那个执行归档的自己。复现的意义从来不在于得到一模一样的数字而在于你亲手把一条因果推断链条完整走过一遍。走完一遍之后你会发现再看任何类似的政策评估论文时那些模型设定、识别假设、稳健性设计都不再是抽象的文字而是一道道你真正操作过的工序。这才是顶刊复现带给你最实在的价值。