连享会   主页 || 推文 || 知乎 || B 站 || 在线课堂

New！ 搜推文，找资料，用 lianxh 命令：
安装： ssc install lianxh, replace
使用： lianxh 合成控制
       lianxh DID + 多期, w

作者： 陈巧雯 (厦门大学); 杜克锐 (厦门大学)
邮箱： chenqiaowen@stu.xmu.edu.cn; kerrydu@xmu.edu.cn

• Title：SFA-xtsfsp：空间随机前沿模型的 Stata 实现
• Keywords：kgitee, xtsfsp, Orea & Álvarez, Galli, 空间随机前沿模型, spatial SFA, stochastic frontier analysis, 效率分析

1. 简介

我们在之前的两篇推文中（杜克锐, 2024, SFA-spsfe：空间随机前沿模型的 Stata 实现；王珞嘉, 2024, SFA-sdsfe：内生空间随机前沿模型的 Stata 实现），介绍了如何使用 spsfe 和sdsfe命令实现空间随机前沿模型估计。spsfe命令实现 Kultu et al. (2020) 提出的内生空间随机前沿模型估计，考虑了无效率项的影响因素和变量内生性。sdsfe命令实现了 Galli (2023)考虑无效率项空间相关性的综合空间随机前沿模型估计，并借鉴了 Kultu et al. (2020)控制函数的思想允许前沿和环境变量内生。作为本系列的第三篇推文，本文将介绍具有多重空间效应的空间随机前沿模型及其 Stata 估计命令xtsfsp，基于 Orea & Álvarez (2019)和 Galli (2022)提出的模型，将无效率项和扰动项的空间相关性以及前沿函数的全局和局部空间溢出纳入随机前沿分析。

2. 模型设定

基于 Wang & Ho（2010）模型的转置，Orea & Álvarez（2019）提出了空间随机前沿模型，该模型考虑了空间相关的无效率项和扰动项。该模型的公式如式(1)-(3)所示，其中 $i=1, \ldots, N；t=1, \ldots, T$：
 

$$Y_{i t}=X_{i t}^{\prime} \beta+\tilde{v}_{i t}-s \tilde{u}_{i t} \tag{1}$$

$$\tilde{v}_{i t}=v_{i t}+\gamma W_i^{v t} \tilde{v}_{. t}\tag{2}$$

$$\tilde{u}_{i t}=u_{i t}+\tau W_i^{u t} \tilde{u}_{. t} \tag{3}$$

 
式(1)描述了随机前沿函数，其中 $Y_{i t}$ 是因变量，$X_{i t}$ 是一个塑造前沿的 $k \times 1$ 变量向量；在生产函数中 $s=1$，在成本函数中 $\mathrm{s}= -1$ 。$\tilde{v}_{i t}$ 和 $\tilde{u}_{i t}$ 分别表示异质性噪音项和无效率项。在式(2)和(3)中，滞后自变量 (即 $\tilde{v}_{. t}$ 和 $\tilde{u}_{. t}$ ) 的空间参数 $W_i^{v t}=$ $\left(W_{i 1}^{v t}, \ldots, W_{i N}^{v t}\right)$ 及 $W_i^{u t}=\left(W_{i 1}^{u t}, \ldots ,W_{i N}^{u t}\right)$ 是两个已知的 $1\times \mathrm{N}$ 截面权重向量，分别刻画了异质性噪声项和无效率项的截面关系结构；$\tilde{v}_{. t}=\left(\tilde{v}_{1 t}, \ldots, \tilde{v}_{N t}\right)^{\prime}$，并且 $\tilde{u}_{. t}=\left(\tilde{u}_{1 t}, \ldots, \tilde{u}_{N t}\right)^{\prime}； v_{i t}$ 为服从分布 $N\left(0, \sigma_{v, i t}^2\right)$ 的随机变量，$u_{i t}=h\left(Z_{i t}^{\prime} \delta\right) u_t^*$ 。 $h\left(Z_{i t}^{\prime} \delta\right)=\sqrt{\exp \left(Z_{i t}^{\prime} \delta\right)}$ 为 scaling 函数，$Z_{i t}$ 为影响个体无效率的变量的 $1 \times 1$ 向量，$u_t^*$ 为服从分布 $N^{+}(0,1)$ 的非负随机变量。与 Orea 和 Álvarez (2019) 的原始设定不同，我们假设 $\sigma_{v, i t}^2=\exp \left(D_{i t}^{\prime} \eta\right)$ 来解释依赖于变量向量 $D_{i t}^{\prime}$ 的异质性误差方差 。更进一步，我们令 $u_t^*$ 的方差等于 1，从而允许 $Z_{i t}$ 包含一个常数。利用矩阵表示法，我们可以重写式(2)和(3)：
 

$$\tilde{v}_{. t}=\left(I_N-\gamma W^{v t}\right)^{-1} v_{. t} \tag{4}$$

$$\tilde{u}_{. t}=\left(I_N-\tau W^{u t}\right)^{-1} h\left(Z_{. t} \delta\right) u_t^*=\tilde{h}_{. t} u_t^* \tag{5}$$

 
其中 $Z_{. t}=\left(Z_{1 t}, \ldots, Z_{N t}\right)^{\prime} ; \tilde{h}_{. t}=\left(I_N-\tau W^{u t}\right)^{-1} h\left(Z_{. t} \delta\right)$ 。

上述模型通过空间自回归（SAR）过程捕捉扰动项和无效率项的空间相关性。

Galli（2022）进一步将因变量和投入变量的空间滞后项纳入 Orea 和 Álvarez（2019）模型，从而测度了影响前沿函数的全局和局部空间溢出。模型表示为：
 

$$Y_{i t}=\rho W_i^{y t} Y_{. t}+X_{i t}^{\prime} \beta+W_i^{x t} X_{. t} \theta+\tilde{v}_{i t}+s \tilde{u}_{i t} \tag{6}$$

 
其中 $W_i^{y t}=\left(W_{i 1}^{y t}, \ldots, W_{i N}^{y t}\right)$ 和 $W_i^{x t}=\left(W_{i 1}^{x t}, \ldots, W_{i N}^{x t}\right)$ 是两个已知的 $1 \times \mathrm{N}$ 截面权重向量 ，分别与空间滞后因变量 $\left(Y_{. t}\right)$ 和自变量 $\left(X_{. t}\right)$ 相关; $Y_{. t}=$ $\left(Y_{1 t}, \ldots, Y_{N t}\right)^{\prime}$ 以及 $X_{. t}=\left(X_{1 t}, \ldots, X_{N t}\right)^{\prime}$ 。

Galli (2022) 对四种不同类型的空间相关性给出了全面的描述：因变量 $Y_{i t}$ 的全局溢出，输入变量 $X_{i t}$ 的局部溢出，异质性噪声 $v_{i t}$ 和无效率 $u_{i t}$ 的横截面相关性。我们称这个完整的模型为 ”yxuv-SAR“。可以对特定的参数施加一些限制，生成以下模型 (见 Table 1)，并通过 xtsfsp命令进行估计。

3. Stata 新命令：xtsfsp

3.1 命令安装

打开 Stata ，在命令窗口输入以下命令并运行，即可完成安装 (请保证安装时网络畅通) 。

net install kgitee, from(https://gitee.com/kerrydu/kgitee/raw/master) replace

kgitee xtsfsp // 安装程序文件，下载范例数据文件到当前工作路径下

若上述指令顺利执行，则 Stata 结果窗口中将显示如下结果：

. kgitee xtsfsp checking xtsfsp consistency and verifying not already installed...
installing into D:\Software\Stata17\ado\plus...
installation complete.
checking xtsfsp consistency and verifying not already installed...

copying into current directory...
      copying xtsfsp\_ex1.dta
      copying xtsfsp\_ex2.dta
      copying xtsfsp\_ex3.dta
      copying xtsfsp\_w1.mmat
      copying xtsfsp\_w2.mmat
      copying xtsfsp\_w3.mmat
ancillary files successfully copied.

这些存储于当前工作路径下的 .dta 和 .mmat 文件将在第 4 小节的 Stata 操作示例中调用。

注意：xtsfsp 依赖于 Jann (2005) 提供的 moremata 软件包。如果尚未安装，可在 Stata 命令窗口输入以下代码进行安装。

ssc install moremata

3.2 命令语法

xtsfsp 的基本语法如下：

xtsfsp depvar [indepvars] [, options]

其中，depvar 为被解释变量，indepvars 为解释变量。options 包括：

• 前沿函数选项

  • noconstant：前沿函数中不包含常数项
  • cost：指定前沿函数为成本函数。默认前沿函数为生产函数
  • wxvars(varlist)：前沿函数中的空间杜宾变量

• 方差函数选项

  • uhet(varlist[, noconstant])：无效率项 scaling 函数的解释变量；使用 noconstant 则函数不包含常数项
  • vhet(varlist[, noconstant])：异质性误差方差函数的解释变量；使用 noconstant 则函数不包含常数项

• 空间权重选项

  • wy(W1 [W2...WT][,mata array])：为空间滞后项指定空间权重矩阵
  • wx(W1 [W2...WT][,mata array])：为前沿函数中的空间杜宾项指定空间权重矩阵
  • wu(W1 [W2...WT][,mata array])：为技术无效率项指定空间权重矩阵
  • wv(W1 [W2...WT][,mata array])：为扰动项指定空间权重矩阵
  • normalize(string)：指定空间权重矩阵的归一化方法

• 回归选项

  • initial(matname)：指定极大似然估计初始值矩
  • endvars(varlist)：指定内生变量
  • iv(varlist)：指定工具变量
  • leaveout(varlist)：在内生变量回归中剔除的变量
  • mlsearch(search_options)：指定搜索初始值的命令
  • delve：基于特定模型进行初始值窥探
  • mlplot：使用 ml plot 找到更好的初始值
  • mlmodel(model_options)：控制 ml 模型的选项
  • mlmax(maximize_options)：控制 ml 最大化的选项
  • delmissing：将缺失值对应的空间权重矩阵元素删除

• 结果报告选项

  • nolog：省略极大似然估计迭代日志的显示
  • genwvars：生成空间杜宾和空间滞后项

4. Stata 操作示例

xtsfsp 命令具有极大的灵活性，可以对各种形式的空间相关性进行各种估计。具体来说，xtsfsp可以指定不同的空间成分组合，还可以选择不同的、随时间变化的空间权重矩阵。此外，xtsfsp还可以指定随机误差的条件异方差。在本节中，我们将介绍两个示例，演示如何使用 xtsfsp命令。所有示例均在在 Stata17, Intel(R) Core(TM) i7-9750H CPU 6 cores @2.60GHz, 8GB RAM @2667MHz, HDD 配置下进行。

4.1 具有不同空间权重矩阵的 xuv-SAR 模型

在这个例子中，我们估计无效率项、扰动项和空间杜宾项具有不同空间权重矩阵的空间随机前沿模型。cost选项指定前沿函数类型为成本函数，uhet(z)选项指定了无效率方差函数中的控制变量为 $z$。矩阵b指定了最大似然估计的初始值。良好的初始值将有利于空间随机模型的拟合，空间相关参数的初始值可以通过使用 frontier 和 sfpanel命令拟合非空间随机模型获得。wu(w2,mata)选项用于指定无效率项的空间权重矩阵，wv(w1,mata)为扰动项指定空间权重矩阵，wxvars(x)指定空间杜宾项的空间权重矩阵。

*-空间权重矩阵读入mata环境
. mata mata matuse xtsfsp_w2, replace
(loading w1[300,300], w10[300,300], w2[300,300],
 w3[300,300], w4[300,300], w5[300,300], w6[300,300],
 w7[300,300], w8[300,300], w9[300,300])


. local w w1 w2 w3 w4 w5 w6 w7 w8 w9 w10

*-导入示例数据
. use xtsfsp_ex2.dta, replace

*-设置初始估计参数
. mat b=(2,0.5,1,-1.5,4,-1.5,0.6,0.6)

. timer clear 1
. timer on 1

*-估计模型
. xtsfsp y x, cost uhet(z) wu(w2,mata)  ///
              wv(w1,mata)  ///
              wxvars(x)    ///
              wx(`w',mata) ///
              init(b) genwvars

warning: genwxvars was replaced by genwvars, see help xtsfsp

initial:       log likelihood = -2206.8046
rescale:       log likelihood = -2206.8046
rescale eq:    log likelihood = -2175.8736
Iteration 0:   log likelihood = -2175.8736
... ...
Iteration 20:  log likelihood =  -1911.966

initial:       log likelihood = -2204.3268
rescale:       log likelihood = -2204.3268
rescale eq:    log likelihood = -2173.5485
Iteration 0:   log likelihood = -2173.5485
... ...
Iteration 26:  log likelihood = -1911.966

Spatial frontier model(xuv-SAR)    Number of obs =    3,000
                                  Wald chi2(2)   = 57442.57
Log likelihood = -1911.966        Prob > chi2    =   0.0000

-----------------------------------------------------------
        y |  Coeff     SE        z   P>|z|       [95% CI]
----------+------------------------------------------------
 frontier |
        x |  1.996   .008   238.00   0.000     1.979  2.012
      W_x |  .5069   .022    22.90   0.000      .464   .550
    _cons |  .9917   .013    77.91   0.000      .967  1.017
----------+------------------------------------------------
 /lnsigv2 | -1.615   .026   -62.04   0.000    -1.666 -1.564
----------+------------------------------------------------
     uhet |
        z |  3.999   .002  1855.94   0.000      3.995 4.004
    _cons | -1.699   .447    -3.80   0.000     -2.576 -.822
----------+------------------------------------------------
       Wv |
    _cons |  .5878   .058    10.06   0.000      .473   .702
----------+------------------------------------------------
       Wu |
    _cons |  .6214   .001   728.34   0.000      .620   .623
----------+------------------------------------------------
      tau |  .3010   .000   776.13   0.000      .3003  .302
    gamma |  .2857   .027    10.65   0.000      .2323  .337
-----------------------------------------------------------

Note: Wv:_cons and Wu:_cons are the transformed parameters;
      gamma and tau are their origin metrics in spatial components, respectively.
      W_(x) represent Spatial Durbin terms W(x)

. timer off 1
. timer list
   1:    178.26 /        1 =     178.2620

为了演示如何使用 delmissing 选项，我们将 $y_{it}$ 的两个观测值替换为缺失值，并增加 delmissing 选项估计上述模型。生成的变量 $e\_sample$ 记录了回归样本。本例的执行时间相对较长，超过 17 分钟，是不包含缺失值模型估计的 12.2 倍。这是因为时变维度以及缺失值导致的数据不平衡，需要计算随时间变化的空间权重矩阵，计算所需时间增加。

*-空间权重矩阵读入mata环境
. mata mata matuse xtsfsp_w2, replace
(loading w1[300,300], w10[300,300], w2[300,300], w3[300,300], w4[300,300], w5[300,300], w6[300,300], w7[300,300], w8[300,300], w9[300,300])

. local w w1 w2 w3 w4 w5 w6 w7 w8 w9 w10

*-导入示例数据
. use xtsfsp_ex2.dta, replace

*-设置初始估计参数
. mat b=(2,0.5,1,-1.5,4,-1.5,0.6,0.6)

*-生成缺失值
. replace y = . if _n==1 | _n==100
(2 real changes made, 2 to missing)

. timer on 2

*-估计模型
. xtsfsp y x, cost uhet(z) wu(w2,mata)  ///
              wv(w1,mata)  ///
              wxvars(x)    ///
              wx(`w',mata) ///
              init(b)      ///
              delmissing

warning: genwxvars was replaced by genwvars, see help xtsfsp
missing values found. The corresponding units are deleted from the spmatrix

initial:       log likelihood = -2204.3268
rescale:       log likelihood = -2204.3268
rescale eq:    log likelihood = -2173.5485
Iteration 0:   log likelihood = -2173.5485
... ...
Iteration 26:  log likelihood = -1909.3843

Spatial frontier model(xuv-SAR)    Number of obs =    2,998
                                   Wald chi2(2)  = 57432.37
Log likelihood = -1909.3843        Prob > chi2   =   0.0000

-----------------------------------------------------------
        y |  Coeff     SE        z   P>|z|       [95% CI]
----------+------------------------------------------------
 frontier |
        x |  1.996   .008   238.01   0.000      1.98  2.012
      W_x |  .5065   .022    22.89   0.000      .463  .5499
    _cons |  .9914   .013    77.97   0.000      .967  1.016
----------+------------------------------------------------
 /lnsigv2 | -1.616   .026   -62.05   0.000    -1.667 -1.564
----------+------------------------------------------------
     uhet |
        z |  3.999   .002  1856.92   0.000     3.995  4.004
    _cons | -1.699   .447    -3.80   0.000    -2.575 -.8223
----------+------------------------------------------------
       Wv |
    _cons |  .5865   .058    10.03   0.000     .472   .7011
----------+------------------------------------------------
       Wu |
    _cons |  .6214   .001   728.88   0.000     .620   .6231
----------+------------------------------------------------
      tau |  .3010   .000   776.71   0.000     .3003  .3018
    gamma |  .2851   .027    10.62   0.000     .2317  .3368
-----------------------------------------------------------

Note: Wv:_cons and Wu:_cons are the transformed parameters;
      gamma and tau are their origin metrics in spatial components, respectively.
      W_(x) represent Spatial Durbin terms W(x)
      Missing values found
      The regression sample recorded by variable __e_sample__

. timer off 2

. timer list
   1:    178.26 /        1 =     178.2620
   2:   2180.13 /        1 =    2180.1270

4.2 具有条件异方差随机误差项的 uv-SAR 模型

本例估算一个具有时变空间权重矩阵和条件异方差随机误差的限制模型 uv-SAR。本例计算时间接近 15 分钟，主要原因是空间权重矩阵的配置随时间变化导致运算较为复杂。

*-空间权重矩阵读入mata环境
. mata mata matuse xtsfsp_w3,replace
(loading w1[300,300], w10[300,300], w2[300,300], w3[300,300], w4[300,300], w5[300,300], w6[300,300], w7[300,300], w8[300,300], w9[300,300])

. local w w1 w2 w3 w4 w5 w6 w7 w8 w9 w10

*-导入示例数据
. use xtsfsp_ex3.dta, replace
. xtset id t

. timer on 3

*-估计模型
. xtsfsp y x, wu(`w',mata) wv(`w', mata) ///
              uhet(z) vhet(d)

warning: genwxvars was replaced by genwvars, see help xtsfsp
searching initial values...
initial:       log likelihood = -17810.346
... ...
Iteration 10:  log likelihood = -5803.0516

initial:       log likelihood = -5803.0516
rescale:       log likelihood = -5803.0516
rescale eq:    log likelihood = -5772.5485
Iteration 0:   log likelihood = -5772.5485
... ...
Iteration 26:  log likelihood = -5803.0516

Spatial frontier model(uv-SAR)    Number of obs =   3,000
                                  Wald chi2(1)  = 7472.38
Log likelihood = -5803.0516       Prob > chi2   =  0.0000

-----------------------------------------------------------
        y |  Coeff     SE        z   P>|z|       [95% CI]
----------+------------------------------------------------
 frontier |
        x |  1.988   .023    86.44   0.000      1.943  2.033
    _cons |   .900   .074    12.20   0.000      .755  1.044
----------+------------------------------------------------
 lnsigv2  |
        d |   .990   .026    38.19   0.000      .939  1.041
    _cons |   .995   .026    38.34   0.000      .944  1.046
----------+------------------------------------------------
     uhet |
        z |  1.019   .046    22.22   0.000      .929  1.109
    _cons |  1.123   .463     2.42   0.015      .215  2.030
----------+------------------------------------------------
       Wv |
    _cons |   .617   .043    14.29   0.000      .532   .702
----------+------------------------------------------------
       Wu |
    _cons |   .647   .070     9.20   0.000      .509   .785
----------+------------------------------------------------
      tau |   .299   .020    15.21   0.000      .260   .337
    gamma |   .313   .032     9.85   0.000      .249   .374
-----------------------------------------------------------

Note: Wv:_cons and Wu:_cons are the transformed parameters
      gamma and tau are their origin metrics in spatial components, respectively.


. timer off 3
. timer list
   1:    178.26 /        1 =     178.2620
   2:   2180.13 /        1 =    2180.1270
   3:    899.26 /        1 =     899.2630

5. 结语

关于xtsfsp命令更详细的描述参见下面这篇工作论文，欢迎读者查阅：

Kerui Du, Luis Orea, Inmaculada C. Alvarez. Fitting spatial stochastic frontier models in Stata

xtsfsp还处于不断完善过程中，命令将不断更新。如使用过程遇到 bugs，请发送邮件至 kerrydu@xmu.edu.cn，感谢大家的宝贵意见。

6. 参考资料

• Galli, F. 2023. A spatial stochastic frontier model introducing inefficiency spillovers. Journal of the Royal Statistical Society Series C: Applied Statistics, 72(2), 346–367. -Link- , -PDF- , -Google-.

• Galli, F., 2022. A spatial stochastic frontier model including both frontier and error-based spatial cross-sectional dependence. Spatial Economic Analysis 18, 239–258. -Link- , -PDF-, -Google-.

• Kutlu, Levent, Kien C. Tran, and Mike G. Tsionas. 2020. A Spatial Stochastic Frontier Model with Endogenous Frontier and Environmental Variables. European Journal of Operational Research 286(1): 389–99. -Link- , -PDF- , -Google-.

7. 相关推文

Note：产生如下推文列表的 Stata 命令为：

lianxh 效率 空间随机前沿

安装最新版 lianxh 命令：  

ssc install lianxh, replace

• 彭甲超, 2021, Stata 效率分析：Simar-Wilson 两阶段半参数 DEA-T303, 连享会 No.703.
• 杜克锐, 2024, SFA-spsfe：空间随机前沿模型的 Stata 实现, 连享会 No.1389.
• 连享会, 2021, 连享会：效率分析-前沿模型, 连享会 No.734.
• 黄晨晨, 2024, Lasso-SFA：基于自适应 Lasso 的面板固定效应随机边界分析, 连享会 No.1374.

资源共享

• 连享会资料 ……
• 在线视频：lianxh-class.cn
  • Stata 33 讲，100 万+ 播放，Stata 入门必备，公开课
  • 直击面板数据模型，10 万+ 播放，白话面板模型，公开课
  • … more …
• 论文复现和数据
  • 学术论文重现网站大全
  • 常用数据库 | 人文社科开放数据库
• 主题分类
  • 所有分类 | Stata教程 | 结果输出 | 绘图 | 数据处理
  • 回归分析 | 面板数据 | 交乘项 | 因果推断 | DID | RDD
  • Markdown | 工具软件 | 机器学习 | 其它
• 热门推文
  • Stata：DID 入门教程
  • 稳健性检验！ | 安慰剂检验！ | 内生性！
  • 面板数据模型 |  聚类调整后的标准误

尊敬的老师 / 亲爱的同学们：
连享会致力于不断优化和丰富课程内容，以确保每位学员都能获得最有价值的学习体验。为了更精准地满足您的学习需求，我们诚挚地邀请您参与到我们的课程规划中来。
请您在下面的问卷中，分享您 感兴趣的学习主题或您希望深入了解的知识领域 。您的每一条建议都是我们宝贵的资源，将直接影响到我们课程的改进和创新。
我们期待您的反馈，因为您的参与和支持是我们不断前进的动力。感谢您抽出宝贵时间，与我们共同塑造更加精彩的学习旅程！https://www.wjx.cn/vm/YgPfdsJ.aspx# 再次感谢大家宝贵的意见！

关于我们

• Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。more……
• 扫码加入连享会微信群，提问交流更方便