Stata：Simar-Wilson两阶段半参数DEA-simarwilson

发布时间：2021-11-22 阅读 2518

Stata连享会主页 || 视频 || 推文 || 知乎 || Bilibili 站

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装：
. ssc install lianxh
详情参见帮助文件 (有惊喜)：
. help lianxh
连享会新命令：cnssc, ihelp, rdbalance, gitee, installpkg

✌ 课程详情： https://gitee.com/lianxh/Course

⛳ 课程主页： https://gitee.com/lianxh/Course

⛳ Stata 系列推文：

全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
计量专题 | 论文写作 | 数据分享 | 专题课程
结果输出 | Stata绘图 | 数据处理 | Stata程序
回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
文本分析-爬虫 | Python-R-Matlab | 机器学习
Markdown | 工具软件 | 其它

☝ PDF下载 - 推文合集

作者：彭甲超 (中国地质大学)
邮箱：pengjiachao@cug.edu.cn

编者按：本文主要摘译自下文，特此致谢！
Source：Badunenko O, Tauchmann H. Simar and Wilson two-stage efficiency analysis for Stata[J]. The Stata Journal, 2019, 19(4): 950-988. -PDF-

1. 简介

DEA 方法作为一种非参数效率评价方法，其生产前沿面和效率值估计依赖于样本数据，因而对样本具有较强的敏感性 ( Simar 和 Wilson，2007)。在利用效率值进行回归时，由于存在序列相关性、影响因素内生性等问题，传统 DEA 两阶段回归的统计推断会失效 (Simar 和 Wilson 2007)。例如，首先使用 DEAP、MATLAB 或 Stata 计算出决策单元效率值，然后使用 tobit 或 reg 等方法进行影响因素回归。

传统割裂的 DEA 两阶段法缺乏清晰的基础数据生成过程理论，且忽略估计的 DEA 效率得分是根据共同的数据样本计算得出，将它们视为独立的观测值会因与序列相关而无效推断 (Simar 和 Wilson，2007)。鉴于此，Simar 和 Wilson (2007) 提出的两阶段半参数 DEA 方法。

该方法描述了一个与两阶段估计程序一致的数据生成过程，提出了两个与假定的数据生成过程一致的 bootstrap 参数程序，有效避免决策单元生成效率值与影响因素回归估计的偏误问题。

2. 计算过程

传统的 DEA 两阶段法将可能会影响决策单元效率的因素统称为环境变量。环境变量与投入、产出变量不同，都是决策单元无法直接控制到的外在影响因素 (冯宗宪等，2011)。Simar and Wilson 两阶段效率分析法分两步进行：

第一阶段根据决策单元的投入和产出变量，应用 DEA 方法估计效率值；
第二阶段再将第一阶段求得的效率值对环境变量进行截断回归分析。上

述步骤均在一个程序内完成，不需要额外重新整理数据。Simar and Wilson 两阶段效率分析法与随机前沿分析过程并不存在明显差异 (易明等，2019)，区别仅在于具体的计算方法选择。

具体来看，Simar and Wilson 两阶段效率分析法在第一阶段，需要使用投入产出数据集测算出产出 (投入) 导向的效率值。以产出导向 VRS 模型为例，对于特定的数据点 $(x_{i}, y_{i})$ 有：

{\hat{δ}}_{i} = \hat{δ} (x_{i}, y_{i} | \hat{p r o d u c t i o n}) = m a x {θ > 0 | θ y_{i} \leq Y q, x_{i} \geq X q, q \in R_{+}^{n}}

其中， $i = 1, \dots, n$ 表示决策单元； $\hat{p r o d u c t i o n}$ 表示生产集合； $Y = (y_{1}, \dots, y_{n}), X = (x_{1}, \dots, x_{n})$ 分别表示 $n \times 1$ 维产出和投入向量； $q$ 是由强度变量组成的 $n \times 1$ 的向量； $θ$ 衡量样本点 $(x_{i}, y_{i})$ 与效率前沿之间的径向距离； $1 \leq {\hat{δ}}_{i} < \infty$ ，而 ${\hat{δ}}_{i} - 1$ 表示当投入量不变时第 $i$ 个决策单元的产出可按比例增加的量； $1 / {\hat{δ}}_{i}$ 表示技术效率值。

在获得效率估计值后的第二阶段，对于 ${\hat{δ}}_{i} > 1$ 的 $s$ 个样本点，建立回归模型：

{\hat{δ}}_{i} = z_{i} β + ξ_{i}

应用最大似然估计可计算环境变量 $z_{i}$ 的系数 $\hat{β}$ ，及误差项标准差的估计值 ${\hat{σ}}_{ϵ}$ 。接着，将以下三个步骤循环进行 L 次以获得一个 bootstrap 自助法估计值的集合 $A = {({\hat{β}}^{*}, {\hat{σ}}_{ϵ}^{*})_{b}}_{b = 1}^{L}$ ，一般情况下 $L = 100$ (Simar 和 Wilson，2007)：

第一步，对每一个 $i = 1, \dots, s$ ，从在点 $(1 - z_{i} \hat{β})$ 处左截尾的 $N (0, {\hat{σ}}_{ϵ}^{2})$ 分布随机抽取 $ϵ_{i}$ ；
第二步，再次地对于 $i = 1, \dots, s$ ，计算 $δ_{i}^{*} = z_{i} \hat{β} + ϵ_{i}$ ；
第三步，使用最大似然估计 $δ_{i}^{*}$ 对 $z_{i}$ 的截尾回归，得出估计值 $({\hat{β}}^{*}, {\hat{σ}}_{ϵ}^{*})$ 。

在获得集合 $A$ 后，就可以对 $β$ 中的每个元素及 $σ_{ϵ}$ 使用集合 $A$ 与 $\hat{β}$ 、 ${\hat{σ}}_{ϵ}$ 构建被估计出来的置信区间。具体统计推断过程及 bootstrap 的模拟结果可参考 Badunenko 和 Tauchmann (2019)。

3. 命令介绍

* 命令安装
. findit simarwilson  // 查找最新版
. ssc install simarwilson, replace
. net get st0585_1.pkg, replace

* 命令语法
simarwilson [(outputs = inputs)] [depvar] indepvars [if] [in] [weight] [, options]

其中，主要变量如下：

outputs 代表生产数据集的产出变量，inputs 代表投入变量，数据均为非负变量，并且投入和产出变量数量不得超过决策单元的数量；
depvar 为指定的被解释变量-效率值，该度量将作为因变量进入回归模型。simarwilson 期望 depvar 是一种径向效率，效率范围在 $(0, 1]$ 区间或 $[1, + \infty)$ 区间。如果 depvar 的某些值小于1，而其他值超过1，则 simarwilson 会发出警告，并根据指定 nounit 的方式忽略观察值；
indepvar 表示解释变量。解释变量和被解释变量均不允许使用带有时间序列运算符等，如 L. 或 F.。

options 主要包括：

algorithm(1|2) 指定命令使用算法 1 还是算法 2。为了计算经偏差校正的效率值，算法 2 涉及 bootstrap DEA 的引导程序。算法 2 要求指定 (output=inoput1 input2 input3)。simarwilson 默认值为算法 1；
notwosided 默认命令 simarwilson 应用截断回归模型，无论效率得分在 $(0, 1]$ 区间还是 $[1, + \infty)$ 区间。对于 $(0, 1]$ 内的效率得分，默认使用双向截断回归模型并从双向截断正态分布中抽样。应用 Simar 和 Wilson(2007) 命令时，只考虑 $[1, + \infty)$ 内的效率值，即在没有指定双向的情况下，simarwilson 第二阶段的回归模型在面向产出角度和面向投入的效率之间没有区别，因而不建议在算法 2 中使用 notwoside；
rts(crs|nirs|vrs) 指定在哪种假设下对所考虑的生产过程的规模报酬进行技术效率度量，crs 要求规模报酬不变，nirs 要求规模报酬递增，而 vrs 要求可变规模收益。默认值为 rts(vrs)，rts()一般通过 teradial 起作用。如果使用外部估计的效率值，则指定 rts() 无效；
base(output|input) 指定技术效率的径向度量的方向或基数。产出径向角度一般设定为 base(output)，而投入径向角度则设定为 base(input)。默认值为base(output)；
invert 的设定代表是否使用 Shephard 代替 Farrell计算技术效率。

4. 实际操作

4.1 数据导入

为了说明 simarwilson 如何使用，作者讨论了国家一级的治理质量 (包括司法系统的质量) 是否对以国内生产总值 (GDP) 为产出变量的技术效率产生影响。该分析基于第 9 版 Penn World Table 数据库和世界经济论坛提供的全球竞争力报告跨国面板数据 (2018-02-26 版)，两个数据库均可以在相关网站上公开获得。

作者另外提供了获取全球竞争力指数 (GCI) 数据的 Stata 命令 gciget，并将导入的数据与 Penn World Table 数据库进行合并。具体来看：

第一步：安装 gciget 命令，下载 GCI Excel 数据并导入。若导入数据失败则需详细参考 Stata 报错的内容并做针对性修改。

. ssc install gciget, replace
. *gciget EOSQ048 EOSQ051 EOSQ144, clear
. gciget EOSQ048 EOSQ051 EOSQ144, clear nowarnings ///
>    url(https://file.lianxh.cn/data/g/gci_dataset_2007-2017.xlsx)

Downloading the gci_dataset_2007-2017.xlsx file
Importing the gci_dataset_2007-2017.xlsx file
(6 vars, 152 obs)

第二步：合并数据，生成 2007-2017 年 131 个国家的跨国面板数据样本。

. *quietly merge 1:1 countrycode year using "https://www.rug.nl/ggdc/docs/pwt90.dta"
. quietly merge 1:1 countrycode year using "https://file.lianxh.cn/data/p/pwt90.dta"

第三步：生成范例数据并使用 teradial 命令快速估计决策单元的 DEA 效率。本文不再对 teradial 命令做介绍，感兴趣的读者可参考 Badunenko 和 Mozharovskyi (2016) 关于 teradial 命令的详细解释。

. net install st0444.pkg, replace
. quietly generate regu = EOSQ048[_n-1] if countrycode == countrycode[_n-1]
. quietly generate prop = EOSQ051[_n-1] if countrycode == countrycode[_n-1]
. quietly generate judi = EOSQ144[_n-1] if countrycode == countrycode[_n-1]
. quietly generate lpop = ln(pop[_n-1]) if countrycode == countrycode[_n-1]
. global g_list "regu prop judi"
. global z_list "regu prop judi lpop c.regu#c.lpop c.prop#c.lpop c.judi#c.lpop"
. set level 90
. set seed 341566575
. teradial rgdpo = ck emp hc if year == 2014 & regu <. & prop <. & judi <. ///
>     & lpop <., tename(te_vrs_o) rts(vrs) base(output) noprint

第四步：查看主要指标的描述性统计结果。

. sum te_vrs_o regu prop judi lpop if e(sample)

    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
    te_vrs_o |        131    1.699949    .6236905          1   5.513838
        regu |        131    3.435143    .6711715   1.846199    5.42263
        prop |        131    4.304648    1.030568   1.610298   6.378975
        judi |        131    3.897085    1.315987   1.113236   6.678279
        lpop |        131    2.566502    1.586448  -1.264066   7.217087

其中，te_car_o 为决策单元的 DEA 效率值，使用 teradial 程序计算得出，效率值范围为 $[1, + \infty)$ ；regu、prop、judi、lpop 分别为参与第二阶段模型回归的环境变量。

4.2 tobit 估计

在第二阶段回归部分，作者使用三种模型分别回归并作对比，这三种模型均适用截断数据类型。此外，在回归结束作者也计算了相关变量的边际效应，以便对比主要变量的计算结果及其 95% 置信区间。

作者首先利用 tobit 模型估算回归系数，tobit 模型的输出表明，有 18个国家的技术效率是完全有效的，而有 113 个国家的技术效率是无效的。估计结果表明司法独立质量并未对技术效率产生显著影响。

. tobit te_vrs_o $z_list, ll(1) nolstretch vsquish

Tobit regression                                    Number of obs     =    131
                                                           Uncensored =    113
Limits: Lower =    1                                    Left-censored =     18
        Upper = +inf                                   Right-censored =      0
                                                    LR chi2(7)        =  20.43
                                                    Prob > chi2       = 0.0047
Log likelihood = -128.7102                          Pseudo R2         = 0.0735
------------------------------------------------------------------------------
    te_vrs_o | Coefficient  Std. err.      t    P>|t|     [90% conf. interval]
-------------+----------------------------------------------------------------
        regu |     -0.393      0.201    -1.95   0.054       -0.726      -0.059
        prop |     -0.520      0.257    -2.02   0.046       -0.947      -0.093
        judi |      0.249      0.189     1.32   0.190       -0.064       0.562
        lpop |     -0.821      0.267    -3.08   0.003       -1.263      -0.379
      c.regu#|
      c.lpop |      0.148      0.069     2.16   0.033        0.035       0.262
      c.prop#|
      c.lpop |      0.125      0.087     1.44   0.153       -0.019       0.270
      c.judi#|
      c.lpop |     -0.086      0.069    -1.24   0.218       -0.201       0.029
       _cons |      4.590      0.775     5.92   0.000        3.306       5.874
-------------+----------------------------------------------------------------
var(e.te_v~o)|      0.410      0.056                         0.327       0.514
------------------------------------------------------------------------------

. quietly margins, dydx($g_list) predict(ystar(1,.)) post
. estimates store tobit

4.3 truncreg估计

作者使用 truncreg 替换 tobit 进行回归，此方法会删除 te_vrs_o = 1 的观测值。估计结果的系数与 tobit 得出的系数相比存在较大差异，但是主要指标符号与 tobit 的估计结果相似。根据 truncreg 的结果，司法质量变量似乎对效率又存在重要影响，因为 judi 及其与 lpop 的交互作用在 10% 的水平上具有统计学意义。但根据 Simar 和 Wilson (2007) 的观点，该结果可能是错误估计的标准误差产物。

. truncreg te_vrs_o $z_list, ll(1) nolstretch vsquish

Truncated regression
Limit: Lower =    1                             Number of obs     =        113
       Upper = +inf                             Wald chi2(7)      =      18.90
Log likelihood = -67.606307                     Prob > chi2       =     0.0085
------------------------------------------------------------------------------
    te_vrs_o | Coefficient  Std. err.      z    P>|z|     [90% conf. interval]
-------------+----------------------------------------------------------------
        regu |     -0.926      0.430    -2.15   0.031       -1.633      -0.219
        prop |     -1.244      0.499    -2.49   0.013       -2.065      -0.423
        judi |      0.778      0.378     2.06   0.039        0.157       1.400
        lpop |     -1.740      0.595    -2.92   0.003       -2.719      -0.761
      c.regu#|
      c.lpop |      0.425      0.172     2.47   0.013        0.142       0.708
      c.prop#|
      c.lpop |      0.258      0.179     1.44   0.150       -0.037       0.553
      c.judi#|
      c.lpop |     -0.259      0.150    -1.73   0.083       -0.506      -0.013
       _cons |      7.448      1.630     4.57   0.000        4.767      10.129
-------------+----------------------------------------------------------------
      /sigma |      0.722      0.093     7.81   0.000        0.570       0.874
------------------------------------------------------------------------------

. quietly margins, dydx($g_list) predict(e(1,.)) post
. estimates store truncreg

4.4 simarwilson 估计

对比 tobit 和 truncreg 的结果发现部分指标显著性存在较大差异，相关统计结论存在矛盾。作者进而使用simarwilson, algorithm(1)。由于外部估计的效率值已可得，因此并未选择 algorithm(2) 重新测算决策单元DEA。

. simarwilson te_vrs_o $z_list, reps(2000)

Simar & Wilson (2007) eff. analysis      Number of obs             =       113
(algorithm #1)                           Number of efficient DMUs  =        18
                                         Number of bootstr. reps   =      2000
                                         Wald chi2(7)              =     21.63
inefficient if te_vrs_o > 1              Prob > chi2(7)            =    0.0029
------------------------------------------------------------------------------
Data Envelopment Analysis:                         externally estimated scores

------------------------------------------------------------------------------
             |  Observed   Bootstrap                           Percentile
inefficiency |      Coef.   Std. Err.      z    P>|z|     [90% Conf. Interval]
-------------+----------------------------------------------------------------
te_vrs_o     |
        regu |  -.9258069   .4021808    -2.30   0.021    -1.617009   -.2838173
        prop |  -1.243902   .4715584    -2.64   0.008    -2.042034   -.5066595
        judi |   .7784162    .356048     2.19   0.029     .1985666     1.37162
        lpop |  -1.739993   .5688841    -3.06   0.002    -2.670876   -.8296476
c.regu#c.l~p |   .4253728   .1611459     2.64   0.008     .1649802    .6971241
c.prop#c.l~p |   .2581352   .1692766     1.52   0.127     -.013744    .5414581
c.judi#c.l~p |  -.2592945   .1400455    -1.85   0.064    -.4829779   -.0310317
       _cons |   7.447817   1.557957     4.78   0.000     4.988534    9.974976
-------------+----------------------------------------------------------------
      /sigma |   .7222912   .0877174     8.23   0.000     .5537709    .8368938
------------------------------------------------------------------------------

. quietly margins, dydx($g_list) post
. estimates store alg_1

同样，作者也使用了 algorithm(2) 进行对比：通过设定产出变量与投入变量命令 (rgdpo = ck emp hc)，以及规模报酬设定 rts(vrs) base(output)，计算偏差校正的效率值，结果保存为 tebc(tebc_vrs_o)。作者选择在偏差校正程序中进行 1000 次抽样，高于 Simar 和 Wilson (2007) 中建议的默认值 algorithm(1)。在此过程中，使用偏差校正后的效率值仅会对估计系数和相关的估计置信区间产生中等程度的影响。

. simarwilson (rgdpo = ck emp hc) $z_list if year == 2014, algorithm(2) rts(vrs) ///
>     base(output) reps(2000) bcreps(1000) tebc(tebc_vrs_o)

Simar & Wilson (2007) eff. analysis      Number of obs             =       131
(algorithm #2)                           Number of efficient DMUs  =         0
                                         Number of bootstr. reps   =      2000
                                         Wald chi2(7)              =     22.67
inefficient if tebc_vrs_o > 1            Prob > chi2(7)            =    0.0019
------------------------------------------------------------------------------
Data Envelopment Analysis:               Number of DMUs            =       131
                                         Number of ref. DMUs       =       131
output oriented (Farrell)                Number of outputs         =         1
variable returns to scale                Number of inputs          =         3
bias corrected efficiency measure        Number of reps (bc)       =      1000
------------------------------------------------------------------------------
             |  Observed   Bootstrap                           Percentile
inefficiency |      Coef.   Std. Err.      z    P>|z|     [90% Conf. Interval]
-------------+----------------------------------------------------------------
tebc_vrs_o   |
        regu |   -.920121    .386252    -2.38   0.017    -1.524188   -.2916807
        prop |  -1.224712   .5001283    -2.45   0.014    -2.044433   -.4245777
        judi |   .6782739   .3637609     1.86   0.062     .0807535    1.257161
        lpop |  -1.811459   .5412532    -3.35   0.001    -2.749651   -.9527789
c.regu#c.l~p |   .4253049   .1469884     2.89   0.004     .1875952    .6568649
c.prop#c.l~p |   .2379338   .1737649     1.37   0.171    -.0451599    .5244849
c.judi#c.l~p |  -.2330191   .1396204    -1.67   0.095    -.4655068   -.0031847
       _cons |   7.937758   1.556489     5.10   0.000     5.416196    10.57548
-------------+----------------------------------------------------------------
      /sigma |   .8780338   .1019137     8.62   0.000     .6827387    1.015325
------------------------------------------------------------------------------

. estimates store alg_2

为了定性地解释结果，作者列示了相关指标估计的平均边际效应。

. estimates table tobit truncreg alg_1 alg_2, title(Estimated mean marginal effects) p

------------------------------------------------------------------
    Variable |   tobit       truncreg      alg_1        alg_2     
-------------+----------------------------------------------------
        regu | -.02001409    .04003719     .0400372    .04087528  
             |     0.8110       0.6720       0.6606       0.6825  
        prop | -.17286701   -.33398801   -.33398801   -.34014625  
             |     0.1211       0.0049       0.0041       0.0123  
        judi |  .02948397    .08804266    .08804266     .0666363  
             |     0.7278       0.3449       0.3318       0.5163  
------------------------------------------------------------------

4.5 特殊情形的估计结果

作者使用 Shephard 距离函数代替 Farrell 距离函数计算技术效率的估计结果。

. simarwilson (rgdpo = ck emp hc) $z_list if year == 2014, algorithm(2) rts(vrs) ///
>     base(output) reps(2000) bcreps(1000) invert

Simar & Wilson (2007) eff. analysis      Number of obs             =       127
(algorithm #2)                           Number of efficient DMUs  =         0
                                         Number of bootstr. reps   =      2000
inefficient if bcscore < 1               Wald chi2(7)              =     89.13
twosided truncation                      Prob > chi2(7)            =    0.0000
------------------------------------------------------------------------------
Data Envelopment Analysis:               Number of DMUs            =       131
                                         Number of ref. DMUs       =       131
output oriented (Shephard)               Number of outputs         =         1
variable returns to scale                Number of inputs          =         3
bias corrected efficiency measure        Number of reps (bc)       =      1000
------------------------------------------------------------------------------
             |  Observed   Bootstrap                           Percentile
  efficiency |      Coef.   Std. Err.      z    P>|z|     [90% Conf. Interval]
-------------+----------------------------------------------------------------
bcscore      |
        regu |   .0662423   .0369831     1.79   0.073     .0037006    .1254749
        prop |   .2044327   .0521238     3.92   0.000     .1223982    .2914231
        judi |  -.1061252    .039003    -2.72   0.007    -.1727934   -.0422705
        lpop |   .2482282   .0503275     4.93   0.000     .1667761    .3319345
c.regu#c.l~p |   -.040327   .0127521    -3.16   0.002    -.0605554   -.0190574
c.prop#c.l~p |  -.0633203   .0172175    -3.68   0.000    -.0914905    -.035442
c.judi#c.l~p |    .052152    .013871     3.76   0.000     .0299023    .0760245
       _cons |  -.3585959   .1468773    -2.44   0.015    -.5991185   -.1174532
-------------+----------------------------------------------------------------
      /sigma |   .1200534   .0078431    15.31   0.000     .1033659    .1293603
------------------------------------------------------------------------------

. quietly margins, dydx($g_list) post
. estimates store alg_2_inv

技术效率值的范围分布于 $(0, 1]$ 和 $[1, + \infty)$ 估计情形。

. simarwilson (rgdpo = ck emp hc) $z_list if year == 2014, algorithm(2) rts(vrs) ///
>     base(output) reps(2000) bcreps(1000) invert notwosided

Simar & Wilson (2007) eff. analysis      Number of obs             =       131
(algorithm #2)                           Number of efficient DMUs  =         0
                                         Number of bootstr. reps   =      2000
inefficient if bcscore < 1               Wald chi2(7)              =     91.41
onesided truncation                      Prob > chi2(7)            =    0.0000
------------------------------------------------------------------------------
Data Envelopment Analysis:               Number of DMUs            =       131
                                         Number of ref. DMUs       =       131
output oriented (Shephard)               Number of outputs         =         1
variable returns to scale                Number of inputs          =         3
bias corrected efficiency measure        Number of reps (bc)       =      1000
------------------------------------------------------------------------------
             |  Observed   Bootstrap                           Percentile
  efficiency |      Coef.   Std. Err.      z    P>|z|     [90% Conf. Interval]
-------------+----------------------------------------------------------------
bcscore      |
        regu |   .0377854    .043478     0.87   0.385    -.0320338    .1098205
        prop |   .2704253   .0561111     4.82   0.000     .1800194    .3617601
        judi |  -.1505534   .0407912    -3.69   0.000    -.2201224   -.0844993
        lpop |   .2511181   .0581711     4.32   0.000     .1544998    .3450201
c.regu#c.l~p |  -.0297073    .014931    -1.99   0.047    -.0559272   -.0061124
c.prop#c.l~p |   -.085326   .0192076    -4.44   0.000    -.1181943   -.0530771
c.judi#c.l~p |   .0694229    .015377     4.51   0.000     .0446474    .0943448
       _cons |  -.4126865   .1685216    -2.45   0.014    -.6948666   -.1438625
-------------+----------------------------------------------------------------
      /sigma |   .1401529   .0085838    16.33   0.000      .120564    .1490527
------------------------------------------------------------------------------

. quietly margins, dydx($g_list) post
. estimates store alg_2_notwo

4.6 图示边际效应

作者的Stata命令也给出了绘制指标的边际效应代码，为节省篇幅不在此处展示，感兴趣的读者可以尝试绘制相关指标的边际效应。

local h_list "$g_list lpop"
foreach h of varlist `h_list'{
   quietly sum `h' if e(sample)
   local mymin = r(min)*0.98
   local myxmin = ceil(`mymin')
   local mymax = r(max)*1.02
   local myxmax = floor(`mymax')
   local mystep = (`mymax'-`mymin')/25
   foreach g of varlist `h_list'{
     local r_list : list h_list - h
     quietly margins if e(sample), dydx(`g') at(`h' = (`mymin' (`mystep') `mymax')   ///
	(asobserved) `r_list')
     quietly marginsplot, xlabel(`myxmin' (1) `myxmax') recast(line) recastci(rarea) ///
	scheme(s2manual)
     quietly graph export "simarwilson2`g'_`h'.eps", as(eps) preview(off) replace    ///
	fontface(Times)
  }
}

5. 结论

simarwilson 实现了Simar 和 Wilson (2007) 的两阶段 DEA 效率分析方法。该命令避免将 DEA 效率值与回归分析过程割裂，使其拥有相对完善的统计基础。该方法通过引入 bootstrap 方法消除传统 DEA 两阶段过程的缺陷，能够有效估算偏差校正的效率值，并通过半参数估计方法引入环境变量参与回归，有别于参数方法估计的随机前沿分析两阶段模型。

6. 参考资料

Simar L, Wilson P W. Estimation and inference in two-stage, semi-parametric models of production processes[J]. Journal of econometrics, 2007, 136(1): 31-64. -PDF-
Simar L, Wilson P W. Two-stage DEA: caveat emptor[J]. Journal of Productivity Analysis, 2011, 36(2): 205-218. -PDF-
Badunenko O, Mozharovskyi P. Nonparametric frontier analysis using Stata[J]. The Stata Journal, 2016, 16(3): 550-589. -PDF- -Slide-
Badunenko O, Tauchmann H. Simar and Wilson two-stage efficiency analysis for Stata[J]. The Stata Journal, 2019, 19(4): 950-988. -PDF-
冯宗宪, 王青, 侯晓辉. 政府投入, 市场化程度与中国工业企业的技术创新效率[J]. 数量经济技术经济研究, 2011, 28(4): 3-17. -Link-
易明, 彭甲超, 吴超. 基于 SFA 方法的中国高新技术产业创新效率研究[J]. 科研管理, 2019, 40(11): 22-31. -Link-

7. 相关推文

Note：产生如下推文列表的 Stata 命令为：
lianxh dea 效率, m
安装最新版 lianxh 命令：
ssc install lianxh, replace

专题：SFA-DEA

专题	嘉宾	直播/回看视频
⭐ 最新专题		文本分析、机器学习、效率专题、生存分析等
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]