一些数据处理的方法

这期推送简单介绍一下我在以往清洗数据的过程中用过的一些野路子。

这期推文其实在上期之后就一直在构思，只是在实际落地的时候有一些小问题需要解决，然后这段时间又在忙其他事情，所以就一直拖到了现在……

$Note:$ 1、下划线字体为链接，可点击跳转；2、推文中的公式与代码块均可左右滑动；3、该文首发于微信公众号DMETP，欢迎关注；4、需要本次推送所使用的数据和代码的朋友，可以在公众号后台对话框内回复关键词coups；5、文中所有外部命令均可通过以下方式获得，以高维固定效应估计命令reghdfe为例：

*- 法一：
   ssc install reghdfe, replace

*- 法二：
   findit reghdfe

*- 法三：
   search reghdfe

< 目录 >

TOC

一、数据格式转换与合并

当我们从EPS中国微观经济数据查询系统按照单年数据查询下载好每一年的csv文件后，假设我们按照年份把这些csv文件分别放到不同的文件夹中，接下来的事情就是把这些csv文件统一转化为dta文件，再将这些同一年份的dta文件纵向append到一个dta文件中。

往期推送的解决方案是：

• 首先，利用批处理对文件重命名；
• 其次，使用StatTransfer软件将csv文件转为dta文件；
• 最后，在Stata中修正乱码并使用for循环进行多个数据集的纵向合并。

整体来看，这样的操作比较繁乱，而且特别容易出错。事实上，Stata有一个csvconvert外部命令能够同时进行多个csv文件的转换与合并，极大地方便了我们的数据清洗工作（除了不能修正乱码）。下面结合工企数据库对该命令进行简要介绍。

首先，在桌面创建一个文件夹，并将其命名为exp，在该文件夹中再创建两个子文件夹：

• 一是raw_data文件夹，raw_data中再分别按照年份创建16个孙文件夹（1998-2013年），分别存放我们下载好的对应年份的原始数据。
• 二是temp_data文件夹，用于存放我们操作过程中产生的缓存数据。

其次，在Stata中定义原始数据及缓存数据存放路径的全局暂元。

global path      C:\Users\KEMOSABE\Desktop\exp
global raw_path  $path\raw_data
global temp_path $path\temp_data

之后，结合循环语句和csvconvert命令将原始数据合并为16个dta文件。

forvalue i = 1998/2013 {
    capture csvconvert $raw_path\\`i' ,                                      ///
            replace output_dir( $temp_path ) output_file(`i')
    }

以1998年为例：

• $raw_path\\1998是1998年所有原始数据存放的路径，在子路径和孙路径中间加两个\的原因是，如果只加一个\，Stata将自动忽略这个符号从而报错。
• replace代表替代该路径下的同名文件。
• output_dir( $temp_path )代表输出的缓存文件的存放路径，这里是C:\Users\KEMOSABE\Desktop\exp\temp_data。
• output_file(1998)代表该缓存文件命名为1998。

整体来看，csvconvert命令确实比之前的方法要简便。但是，一个不能忽视的问题是，在我们得到这16年的dta文件之后，任意打开一个文件可以发现乱码的现象是普遍存在的，包括所有变量名称、字符型数据和标签。之前的方法可以使用Stata自带的转码命令进行转码，但是通过csvconvert输出的数据集却只能对标签进行转码，我尝试了几乎所有能找到的方法，但这些方法都不能奏效。针对这个顽疾，欢迎各位读者在公众号后台与我交流！

二、OP法计算企业TFP

OP法测算TFP的常用外部命令是opreg，在该命令下需要企业退出变量Exit，Exit的取值规则是：若企业在下一年退出市场且于观察期内不复进入时取1，否则取0；特别地，当企业在观测期末年仍旧存续时取0。以工企数据库为例，存在以下四种情况：

• 情况一，企业只有单年观测值（ singleton ），也就是说，某企业在1998-2013年这16年的观测区间内只有一年观测值。对于这种样本，Exit的取值情况不影响回归结果，因为在参与回归时单年观测值将被自动剔除（除非强行不剔除，如reghdfe命令下使用keepsingletons选择项，但这样的后果是统计显著性有偏）。在这里，我们将其取1。
• 情况二，企业存在两年及以上观测值，并且这些观测值在时间上连续，如某企业在2001、2002和2003年这三年内存续，并且在2003年以后不存在（无论其原因是退出市场还是数据本身的缺陷）。对于这种情况，我们令Exit在2001年和2002年取0，2003年取1。
• 情况三，企业存在两年及以上观测值，并且这些观测值在时间上不连续，如某企业在2001、2002和2004年这三年内存续，并且在2004年以后不存在。在这种情况下，我们令Exit在2001年和2002年取0，2004年取1。对于间断的年份（如这里的2003年），我们认为企业在这些间断年份内仍旧保持存续，数据缺失不表明企业不存续。
• 情况四，特别地，如果某企业在2013年存在观测值，由于我们无法得知企业在2014年的存续状态，因此我们令Exit在2013年取0。

我们以一个手工生成的数据集为例。

clear all
input id year emp
      1  2011 1
      2  2008 0
      2  2009 0
      2  2010 0
      2  2011 0
      2  2012 0
      2  2013 0
      3  2008 0
      3  2010 0
      3  2011 0
      3  2012 0
      3  2013 0
      4  2008 0
      4  2011 0
      4  2012 0
      4  2013 0
      5  2008 0
      5  2009 0
      5  2011 0
      5  2012 0
      5  2013 0
      6  2008 0
      6  2009 0
      6  2010 0
      6  2011 0
      6  2012 1
      7  2006 0
      7  2007 1
end
list, sepby(id)

其中，变量emp是企业实际的存续状况，其取值与后面创建的Exit变量保持一致。可以看到，个体1符合上文中的情况一；个体2符合情况二和情况四；个体3、4和5符合情况三和情况四；个体6和个体7符合情况二。

下面是创建Exit变量的代码。

xtset id year
bysort id: gen Exit = (_n == _N)
replace Exit = 0 if year == 2013
bro

可以看到，Exit的取值与emp保持一致。

事实上，opreg命令的编写者Yasar et al.（2008；2012）也提供了变量Exit的生成方法。下面是参照Yasar et al.（2008；2012）的方法生成工企数据库企业退出变量exit的代码。

2 Yasar M, Raciborski R, Poi B. Production Function Estimation in Stata Using the Olley and Pakes MethodJ. Stata Journal, 2008, 8(02): 221-231.

3 Software UpdatesJ. Stata Journal, 2012, 12(01): 167-167.

xtset id year
sort id year
by id: gen count = _N
gen survivor = (count == 16)
gen has98 = 1 if year == 2013
sort id has98
by id: replace has98 = 1 if has98[_n-1]  == 1
replace has98 = 0 if has98 == .
sort id year
by id: gen has_gaps = 1 if year[_n-1] != year-1 & _n != 1
sort id has_gaps
by id: replace has_gaps = 1 if has_gaps[_n-1] == 1
replace has_gaps = 0 if has_gaps == .
sort id year
by id: gen exit = (survivor == 0 & has98 == 0 & has_gaps != 1 & _n == _N)
replace exit = 0 if exit == 1 & year == 2013
drop count survivor has98 has_gaps
bro

可以看到，emp、Exit和exit保持一致。

此外，部分研究者在生成Exit变量时，直接将存续年份不连续的样本剔除，关于这种做法的理论与文献基础我暂时没有找到，因此其合理性存疑。下面的实现代码借鉴了黄河泉老师（2021）在经管之家论坛的回答。若有必要，这段代码在Exit生成之前运行即可。

xtset id year
gen d = 1
tsfill
bysort id: egen tem1 = count(d)
bysort id:  gen tem2 = _N
keep if tem1 == tem2
drop d tem1 tem2 // 剔除存续序列不连续的样本，该步骤的必要性与合理性存疑

OP法测算企业TFP，除了上文介绍的opreg命令，prodest这个外部命令也能实现，且不需要加入Exit变量。prodest除了可以使用OP法测算TFP，还可以使用LP法、WRDG法、MR法、以及OP法和LP法的ACF修正法，具体方法与使用细节可以键入help prodest进行了解。

值得一提的是，在实际测算中，OP、LP、WRDG和MR方法测算的TFP基本保持一致，但ACF法测算的结果与前面几种相比存在较大偏差。篇幅所限，这几种测算方法的结果对比推文没有贴出来，而是放在了网盘中，公众号后台回复关键词coups即可获取下载链接。

三、截取字符串生成新变量

假设有一个字符型变量，在该变量下的数据为字母与数字的混合组合，且长度不一致。现在的需求就是，生成两个新变量，其中一个变量记录字符型变量的字母部分；另外一个变量记录字符型变量的数字部分，且将变量类型改为数值型。

最初的想法是使用substr()函数，但问题是字符型变量的长度不一致，而substr()函数只能截取固定位置的信息，也就是说，substr()函数要求变量的长度保持一致。这里将尝试使用正则表达式（regular expression）。

还是以一个手工生成的数据集为例。

clear all
input str8 x
           "abc"
           "123"
           "abc123"
           "123abc"
           "a1b2c3"
           "1a2b3c"
end

可以看到，该数据集的唯一变量x有多个字符串组合，包括纯字母、纯数字、前字母后数字、前数字后字母以及其他组合类型。下面的实现代码专门针对前三种组合，至于后面几种组合，利用正则表达式也能实现。

gen x1 = regexs(0) if regexm(x, "^([a-z]*)([0-9]*)$")
gen x2 = regexs(1) if regexm(x, "^([a-z]*)([0-9]*)$")
gen x3 = regexs(2) if regexm(x, "^([a-z]*)([0-9]*)$")

destring x3, replace

list, separator(6)
bro

/*
     +-----------------------------+
     |      x       x1    x2    x3 |
     |-----------------------------|
  1. |    abc      abc   abc     . |
  2. |    123      123         123 |
  3. | abc123   abc123   abc   123 |
  4. | 123abc                    . |
  5. | a1b2c3                    . |
  6. | 1a2b3c                    . |
     +-----------------------------+
*/

可以看到，x2和x3就是我们需要的变量。

下面对代码进行解释，由于这里只是正则表达式的简单应用，因此不会做过多分析，但若能够熟练使用该方法，在处理数据时帮助会很大。键入help string_functions详细了解。

• regexs(#)提示截取信息的位置，#等于0说明截取regexm()中逗号前的内容；#等于1说明截取逗号后第一个括号内的内容；#等于2说明截取逗号后第二个括号内的内容。
• regex()提示信息截取的具体内容，双引号内^表示开始，$表示结束；[a-z]表示字母，[0-9]表示数字；*表示截取0个及以上的数字或字母。

四、永续盘存法测算投资额

使用OP法估计企业TFP时，须引入企业投资额作为不可观测生产率冲击的代理变量，以缓解模型中可能存在的同时性偏差（Simultaneity）（Olley & Pakes，1996）。

5 Olley G S, Pakes A. The Dynamics of Productivity in the Telecommunications Equipment IndustryJ. Econometrica, 1996, 64(06): 1263-1297.

因此，在使用OP法估计工业企业TFP时，固定资产投资额是一个必须变量，但在工企数据库中，不存在固定资产投资额字段。一个常见的方法是使用永续盘存法（Perpetual Inventory Method，PIM）间接测算，公式如下：

$$

I{it} = K{it} - (1 - \delta)K_{it-1}

\tag{1}

$$

有两点比较关键的说明：

• 第一，$K{it}$为$i$企业在$t$年的固定资产合计或固定资产净值（这里选择固定资产合计），$(1)$式说明上一年观测值缺失的样本无法计算当年的$I{t}$，因为上一年的固定资产合计数据无法获得，但对于在整个观测区间内至少存在两年观测值的样本（仅有单年观测值的样本在参与回归时将自动被剔除），$t-1$年的企业固定资产合计$K{it-1}$可以以企业固定资产的平均增长率进行估算，假定平均增长率为$g_i$，则$K{it-1}=K{it}/(1+g_i)$（张天华和张少华，2016）。在这里，企业$i$在$t$年的固定资产增长率$g{it}$计算公式如下式$(2)$，其中$t$与$t'$在年份不连续的情况下其差值不等于1，$gi$即为$g{it}$的组内均值。

$$

g{it} = \frac{K{t}-K{t'}}{K{t'}(t-t')}

\tag{2}

$$

• 第二，$\delta$为经济折旧率，对$\delta$的选择是充满有争议的。根据经济学领域的常用选择，可以将经济折旧率设定为9.6%（张军等，2004；何兴强等，2014）。但存在的问题是，张军等（2004）测算的9.6%是各省固定资本形成总额的经济折旧率，也就是说，该数值没有细分到具体地区，也不是我们关注的制造业行业，更关键地，该折旧率的时间跨度为1952-2000年。因此，为了增强估计结果的稳健性，还应该参考其他文献的做法来设定$\delta$，如5%（邵宜航等，2013）、10%（肖文和林高榜，2014；李青原和章尹赛楠，2021）和15%（聂辉华等，2012）。

6 张天华, 张少华. 中国工业企业全要素生产率的稳健估计J. 世界经济, 2016, 39(04): 44-69.

7 张军, 吴桂英, 张吉鹏. 中国省际物质资本存量估算:1952—2000J. 经济研究, 2004(10): 35-44.

8 何兴强, 欧燕, 史卫, 刘阳. FDI技术溢出与中国吸收能力门槛研究J. 世界经济, 2014, 37(10): 52-76.

9 邵宜航, 步晓宁, 张天华. 资源配置扭曲与中国工业全要素生产率——基于工业企业数据库再测算J. 中国工业经济, 2013(12): 39-51.

10 肖文, 林高榜. 政府支持、研发管理与技术创新效率——基于中国工业行业的实证分析J. 管理世界, 2014(04): 71-80.

11 李青原, 章尹赛楠. 金融开放与资源配置效率——来自外资银行进入中国的证据J. 中国工业经济, 2021(05): 95-113.

12 聂辉华, 江艇, 杨汝岱. 中国工业企业数据库的使用现状和潜在问题J. 世界经济, 2012, 35(05): 142-158.

下面以一个手动生成的数据集为例。

clear all

input id year
      1  2011
      2  2008
      2  2009
      2  2010
      2  2011
      2  2012
      2  2013
      3  2008
      3  2011
      3  2012
      3  2013
      4  2008
      4  2013
      5  2008
      5  2009
      5  2011
      5  2012
      5  2013
      6  2008
      6  2009
      6  2010
      6  2011
      6  2012
      7  2006
      7  2007
end
set seed 2021
gen asset = abs(rnormal())

可以看到，个体1只有单年观测值，个体3、4和5的样本年份不连续。如果按照一般的方法使用PIM测算企业固定资产投资额，全部个体在样本初年的指标均无法计算，且年份不连续的个体在年份开始初年的指标也无法计算，如个体3在2011年的固定资产投资额无法计算。

sort id year
bys id: gen growth = (asset[_n] - asset[_n-1]) / (asset[_n-1] * (year[_n] - year[_n-1]))

xtset id year
tsfill
tsappend, add(1)
bys id: replace year = year[1] - 1 if _n == _N
sort id year
bys id: egen growth_mean = mean(growth)

bys id: gen ceiling = 1 if _n == _N
gen asset_hypo = asset
bys id: replace asset_hypo = asset_hypo[_n+1] / (1 + growth_mean)  ///
           if asset == . & asset[_n+1] != . & ceiling != 1

*- 经济折旧率δ = 9.6%（张军等，2004）
bys id: gen investment1 = asset_hypo - asset_hypo[_n-1] * (1 - 0.096)

*- 经济折旧率δ = 5%（邵宜航等，2013）
bys id: gen investment2 = asset_hypo - asset_hypo[_n-1] * (1 - 0.05)

*- 经济折旧率δ = 10%（肖文和林高榜，2014；李青原和章尹赛楠，2021）
bys id: gen investment3 = asset_hypo - asset_hypo[_n-1] * (1 - 0.1)

*- 经济折旧率δ = 15%（聂辉华等，2012）
bys id: gen investment4 = asset_hypo - asset_hypo[_n-1] * (1 - 0.15)

label var investment1 "经济折旧率取9.6%时测算的固定资产投资额"
label var investment2 "经济折旧率取5%时测算的固定资产投资额"
label var investment3 "经济折旧率取10%时测算的固定资产投资额"
label var investment4 "经济折旧率取15%时测算的固定资产投资额"

drop growth growth_mean ceiling asset_hypo
drop if asset == .

list, sepby(id)
*- 只有单年观测值的个体参与回归时将被剔除

sum invest*
#delimit ;
twoway (kdensity investment1)
       (kdensity investment2)
       (kdensity investment3)
       (kdensity investment4)
       ;
#delimit cr

bro

其中，growth为$g{it}$，growth_mean为$g{it}$的组内均值$g_i$，四个investment#即是不同经济折旧率下测算的固定资产投资额。

下面以核密度曲线图来直观呈现这四个investment#的偏离程度，如图 1。可见，这四个investment#的核密度曲线基本重合，也就是说，对于手动生成的数据集，在$5,15$这个区间选择任意一个折旧率进行测算结果均基本保持一致。

图 1 investment#核密度曲线图

五、货币型变量的指数平减

阅读以下文字前，建议把如何对变量进行指数平减弄懂，可参考知乎用户『无宇』个人主页的第一个回答。

这里以一份手工生成的数据集为例，该数据集包含两个个体1997-2019年的工业增加值，这两个个体所属省份均为安徽省（二位数行政区划代码为34），工业增加值以当年价格计算，因此需要进行指数平减，平减工业增加值的常用指数为工业品出厂价格指数PPI。

version 17
clear all

set obs 20
set seed 2021
gen id = 1 in 1/10
replace id = 2 in 11/20
bys id: gen year = 1996 + _n
gen province = 1
gen VA = abs(rnormal())
label var VA "工业增加值(当年价格)"

对工业增加值进行指数平减的大致思路可以分为两步：

• 第一步，将上年=100的PPI转为1997=1；
• 第二步，利用PPI(1997 = 1)对货币型变量进行指数平减。

由于指数平减涉及两个数据集的数据处理与合并，因此下面将使用框架（frame）进行操作。当然，对单个数据集分别进行处理，然后再横向merge到一个数据集也是可行的做法。

需要说明的是，frame是Stata 16新加入的功能，Stata 15及以下版本无法使用。Stata中frame的功能类似于Excel的工作表sheet，方便在同一个操作窗口中打开多份数据集并对数据集进行处理，而不需另外加载Stata软件打开数据集。

先简单介绍一下Stata中的frame系列命令，具体信息请键入help frame进行了解。

• frame dir：显示内存中所有的框架；
• frame pwf：显示当前正在工作的框架，pwf和frame作用与之相同；
• frame create newfraname：创建一个新的框架，并命名为newfraname，mkf newfraname作用与其一致；
• frame drop franame：删除框架franame，但不可删除当前工作框架；
• frame change franame：将工作框架转换到franame，cwf franame作用与之相同；
• frame rename oldfraname newfraname：重命名框架；
• frame franame: command或frame franame {command}在框架franame下运行命令；
• frlink：通过连接变量linkvar建立两个框架之间的连接；
• frget: 通过连接变量将其他框架中的数据复制到当前工作框架中，一般与frlink配合使用，两者的配合使用功能类似于merge命令进行横向合并（append）。

mkf ppi
frame ppi {
    input province year ppi
          34       1997 99.4
          34       1998 96.4
          34       1999 92.9
          34       2000 98.9
          34       2001 98.6
          34       2002 99.82
          34       2003 103.5
          34       2004 108.16
          34       2005 103.25
          34       2006 103.07
          34       2007 103.6
          34       2008 108.4
          34       2009 92.8
          34       2010 108.98
          34       2011 108.3
          34       2012 98.26
          34       2013 98.16
          34       2014 97.4
          34       2015 93.94
          34       2016 98.5
          34       2017 108
          34       2018 103.03
          34       2019 100.3
    end

    replace ppi = ppi / 100
    label var ppi "工业品出厂价格指数(上年 = 1)"
    sort province year

    gen ppi97 = ppi
    bys province: replace ppi97 = 1 if _n
    bys province: replace ppi97 = ppi97[_n-1] * ppi if _n != 1
    label var ppi97 "工业品出厂价格指数(1997 = 1)"

    format ppi* %6.4f
    list
    }

frame dir
pwf

frlink m:1 province year, frame(ppi) gen(linkid)
frget ppi97, from(linkid)

gen VA1 = VA / ppi97
label var VA1 "工业增加值(可比价格)"

bro

其中，linkid是两个框架之间建立联系的关键连接变量，ppi97是以1997年为基期的PPI，VA1是经指数平减后的工业增加值，因此该变量在不同年份可比。