1 数据说明
工业企业数据库每年包含了数十万条企业样本数据，是研究中国工业企业发展的重要微观数据！但是工业企业数据库样本、变量非常多，且不同年份指标还存在差异，给分析带了较大困难！这里我们已2011-2013年3年的工业企业数据库为例来演示工业企业全要素生产率的计算！
参考相关文献，企业投入产出的指标做如下近似处理
-
产出增加值-用主营业务收入近似表示
-
资本-用固定资产合计近似表示
-
总产出-用工业总产值近似表示
-
劳动投入-用就业人口近似表示
-
中间投入-用总产出减去产出增加值近似表示
这里采用Stata16软件来分析。2011-2013年工业企业数据库已经匹配并保存成stata数据格式！（数据及变量中文对照文件放置于文末！）
注意:在计算全要素生产率时，产出数据有两种形式，一种是增加值，一种是总产值。如果用阐述产出用增加值，那么生产函数中不需要加入中间投入变量，但如果产出用总产值，那么生产函数中需要加入中间投入变量！
2 数据预处理
///第一部分：加载工业企业数据.dta文件
use D:\Desktop\工业企业数据.dta
///第二部分：删除缺失值及不符合常理的样本
//2.1 删除工业总产值、销售产值、资产总计、固定资产合计、从业人数为负、小于零或缺失的企业
drop if (indo==0 | indo<0 )
drop if (isov==0 | isov<0 )
drop if (toasA==0 | toasA<0 )
drop if (tofa==0 | tofa<0 )
drop if (fias==0 | fias<0 )
drop if (nvoe==0 | nvoe<0 )
//2.2 删除企业成立时间在1949年之前的企业
drop if opye<1949
//2.3 删除企业销售利润率超过100%的企业
drop if (oppr/isov)>1
//2.4 删除资产总计小于固定资产的企业
drop if toasA<tofa
//2.5 剔除最大和最小的1%变量
egen p01=pctile(nvoe) ,p(1)
g temp=0
replace temp=1 if nvoe<=p01
bys id : egen x=total(temp)
drop if x>=1
drop p01 temp x
egen p99=pctile(nvoe) ,p(99)
g temp=0
replace temp=1 if nvoe>=p99
bys id : egen x=total(temp)
drop if x>=1
drop p99 temp x
///第三部分：计算工业企业投入和产出变量
//3.1 删除id和year重复的样本
duplicates drop id year,force
//3.2 对企业以现价表示的投入和产出数据进行价格平减
g ppi_index=100 if year==2011 /*用ppi指数对工业总产值、工业销售产值进行平减*/
replace ppi_index=98.3 if year==2012
replace ppi_index=96.4 if year==2013
g assets_index=100 if year==2011 /*用固定资产价格指数对固定资产进行平减*/
replace assets_index=101.1 if year==2012
replace assets_index=101.4 if year==2013
g real_add=mabi/ppi_index*100 /*用主营业务收入近似表示产出增加值*/
g real_k=tofa/assets_index*100 /*用固定资产合计近似表示固定资产*/
g real_total=indo/ppi_index*100 /*用工业总产值近似表示总产出*/
//3.3 投入产出变量取对数处理
g lny=ln(real_add)
g lnl=ln(nvoe) /*就业人口近似表示劳动投入*/
g lnk=ln(real_k)
g lnm=ln(real_total-real_add) /*总产出减去产出增加值近似表示中间投入*/
g lnt=ln(real_total)
//3.4 设置面板数据格式
xtset id year
3 OP法计算全要素生产率
计算企业退市变量exit。（数据为2011-2013年共3年，对于其他年份的数据，可以替换其中的"2011"、"2013"和"3"）：
drop if orco==""
xtset id year
gen firmid = id
sort firmid year
by firmid:gen count=_N
gen survivor=count==3
gen has2013=1 if year==2013
sort firmid year
by firmid:replace has2013=1 if has2013[_n-1]==1
replace has2013=0 if has2013==.
sort firmid has2013
by firmid:gen has_gaps=1 if year[_n-1] !=year-1 & _n !=1
sort firmid has_gaps
by firmid:replace has_gaps=1 if has_gaps[_n-1]==1
replace has_gaps = 0 if has_gaps==.
sort firmid year
by firmid:generate exit=survivor==0 & has2013==0 & has_gaps !=1 & _n==_N
replace exit =0 if exit ==1 & year ==2013
生成投资变量：
sort id year
g lni=ln(real_k-0.1*l.real_k)
根据永续盘存法来估算投资，其中折旧设为10%。2011年的样本将全部丢失。
运行opreg命令：
opreg lny, exit(exit) state(lnk) proxy(lni) free(lnl) cvars(year) vce(bootstrap, seed(1) reps(5))
结果会报出如下错误！
该错误主要是因为exit变量中退市企业数据缺失引起的。由于在计算投资变量时，2011年的样本全部缺失，只剩下2012年和2013年的样本，分别运行如下命令：tabstat exit if year==2012,by(exit) s(n)和tabstat exit if year==2013,by(exit) s(n)，结果发现2013年退市企业数为0，因而只有2012年存在退市企业数据（exit=1）。
为了能够正常运行opreg命令，我们对2014年可能退市的企业进行预测！
2013年利润总额小于0的企业，预测其2014年将退市：
replace exit=1 if totp<0 & year==2013
再次运行opreg命令，估计增加值形式下的企业全要素生产率：
opreg lny, exit(exit) state(lnk) proxy(lni) free(lnl) cvars(year) vce(bootstrap, seed(1) reps(5))
predict tfp_op1
该命令采用了产出为增加值的形式，其中exit(),state(),proxy(),free()为必填字段；cvars是op方法两阶段中的共同控制变量，可以不加！
再次运行opreg命令，估计总产值形式下的企业全要素生产率：
opreg lnt, exit(exit) state(lnk) proxy(lni) free(lnl lnm) cvars(year) vce(bootstrap, seed(1) reps(5))
predict tfp_op2
对于总产值形式，生产函数中需要加入中间投入变量（lnm）