Stata: 空间计量模型溢出效应的动态呈现

原文：How to create animated graphics to illustrate spatial spillover effects

中文版PPT 在线浏览

6 March 2018

作者：Di Liu, Senior Econometrician

Stata连享会精品专题 || 精彩推文

这篇文章展示了如何创建动画图形，以展示空间自回归（SAR）模型生成的空间溢出效应。阅读本文后，您可以创建如下的动态图形。

后文包括三个部分：首先，介绍如何估计 SAR 模型的参数；其次，解释 SAR 模型为什么可以产生空间溢出效应；最后，展示如何创建一个动态图形来说明空间溢出效应。

SAR 模型

我想分析一下德克萨斯州各县的凶杀率与失业率的关系。我猜想一个县的凶杀率会对邻近县的凶杀率造成影响。

我想回答两个问题：

如何构建这样一个模型，在这个模型中明确的允许一个县的凶杀率依赖于邻近县的凶杀率。
根据这个模型，如果“达拉斯”（德克萨斯州的一个城市）的失业率上升到 10％，那么与达拉斯县邻近的县的凶杀率将如何变化呢？

建立SAR模型

首先，考虑标准线性模型，即某一个县 iii 的凶杀率 (hratei)({\bf{hrate}}_{i})(hratei)是该县失业率(unemploymenti)({\bf unemployment}_{i})(unemploymenti)线性函数

hratei=β0+β1unemploymenti+ϵi{\bf hrate}_i = \beta_0 + \beta_1 {\bf unemployment}_{i} + \epsilon_i hratei=β0+β1unemploymenti+ϵi

SAR 模型构建了某一个县 iii 的凶杀率 (hratei)({\bf hrate}_i)(hratei) 依赖于其邻近县的凶杀率。此时，需要一些新的符号来定义 SAR 模型。具体的，如果区域 jjj 和区域 iii 相邻，则设 Wi,jW_{i,j}Wi,j 为正数，如果两者不相邻，则 Wi,jW_{i,j}Wi,j 为 0 。由于区域 iii 不可能和自己相邻，所以当 i=ji = ji=j 时，Wi,jW_{i,j}Wi,j 也为 0。

有了这个符号，构建某一个县 iii 的凶杀率依赖于其邻近县凶杀率的 SAR 模型可以写成

hratei=γ1∑j=1NWi,jhratej+β1unemploymenti+β0+ϵi{\bf hrate}_i = \gamma_1\sum_{j=1}^N W_{i,j} {\bf hrate}_{j} + \beta_1 {\bf unemployment}_{i} + \beta_0 + \epsilon_i hratei=γ1j=1∑NWi,jhratej+β1unemploymenti+β0+ϵi

其中，(Wi,j)(W_{i,j})(Wi,j) 定义了第 iii 个县和第 jjj 个县的邻近程度。∑j=1NWi,jhratej\sum_{j=1}^N W_{i,j} {\bf hrate}_{j}∑j=1NWi,jhratej 是县 iii 的邻近县的凶杀率的加权总和，它表示邻近县的凶杀率对县 iii 凶杀率的影响。

将每个县 iii 的邻近县的信息按 (Wi,j)(W_{i,j})(Wi,j) 叠加得到一个矩阵 W{\bf W}W，矩阵W{\bf W}W记录了每个县iii的邻近县的信息，矩阵W{\bf W}W称为空间加权矩阵。

我们所使用的空间加权矩阵具有特殊的结构;每个元素要么是值ccc ，要么是 0，其中 c>0c > 0c>0。这种空间加权矩阵称为 归一化邻近矩阵。

在 Stata 中，我们使用 spmatrix 命令来创建空间加权矩阵，并使用 spregress 命令来拟合截面 SAR 模型。

首先，我从 Stata 网站下载美国各县谋杀率的数据集，并创建一个仅保留德克萨斯州各县数据的子样本。

. /* Get data for Texas counties' homicide rate */
. copy http://www.stata-press.com/data/r15/homicide1990.dta ., replace. use homicide1990
(S.Messner et al.(2000), U.S southern county homicide rates in 1990). keep if sname == "Texas"
(1,158 observations deleted). save texas, replace
file texas.dta saved

直观地，指定所有地域的地图边界的文件被称为shape文件。我们需要从 Stata 网站下载与数据集texas.dta 相匹配的shape数据集 homicide1990_shp.dta，这个数据集记载了德克萨斯州所有县的地图边界信息。然后，我们使用spset命令，将数据集 texas.dta 与 shape 数据集 homicide1990_shp.dta 相关联，即将其设定为空间数据。

. /* Get data for Texas counties' homicide rate */
. copy http://www.stata-press.com/data/r15/homicide1990_shp.dta, replace. spsetSp dataset texas.dtadata:  cross sectionalspatial-unit id:  _IDcoordinates:  _CX, _CY (planar)linked shapefile:  homicide1990_shp.dta

然后使用 spmatrix 命令创建一个标准化的空间权重矩阵。

. /* Create a spatial contiguity matrix */
. spmatrix create contiguity W

根据上述数据和空间权重矩阵，可以估计模型参数。

. /* Estimate SAR model parameters */
. spregress hrate unemployment, dvarlag(W) gs2sls(254 observations)(254 observations (places) used)(weighting matrix defines 254 places)Spatial autoregressive model                    Number of obs     =        254
GS2SLS estimates                                Wald chi2(2)      =      14.23Prob > chi2       =     0.0008Pseudo R2         =     0.0424------------------------------------------------------------------------------hrate |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
hrate        |
unemployment |   .4584241    .152503     3.01   0.003     .1595237    .7573245_cons |   2.720913   1.653105     1.65   0.100    -.5191143    5.960939
-------------+----------------------------------------------------------------
W            |hrate |   .3414964   .1914865     1.78   0.075    -.0338103    .7168031
------------------------------------------------------------------------------
Wald test of spatial terms:          chi2(1) = 3.18       Prob > chi2 = 0.0745

空间溢出

现在我们准备解决第二个问题，利用上文提到的 spregress 命令估计的结果计算空间溢出效应，具体可以分三步进行：

利用原始数据预测谋杀率。
将达拉斯（“Dallas”）的失业率改为 10％并再次预测谋杀率。
计算两次预测之间的差异并绘图。

. preserve /* save data temporarily */. /* Step 1: predict homicide rate using original data */
. predict y0
(option rform assumed; reduced-form mean). /* Step 2: change Dallas unemployment rate to 10%, and predict again*/
. replace unemployment = 10 if cname == "Dallas"
(1 real change made)
. predict y1
(option rform assumed; reduced-form mean). /* Step 3: Compute the prediction difference and map it*/
. generate double y_diff = y1 - y0
. grmap y_diff, title("Global spillover")
. restore /* return to original data */

上面的图表显示，达拉斯失业率的变化，不仅影响了达拉斯本地区的谋杀率，也影响了与达拉斯相邻的地区的谋杀率。也就是说，达拉斯的变化波及到附近的县，这种影响被称为溢出效应。

SAR 模型与空间溢出

在本节中，我将说明为什么SAR模型会产生溢出效应。在此过程中，我提供了一个用于创建动画图形的效果公式。
SAR 模型的矩阵形式是

y=λWy+Xβ+ϵ{\bf y} = \lambda {\bf W} {\bf y} + {\bf X}\beta + \epsilon y=λWy+Xβ+ϵ

求解y{\bf y}y的值

y=(I–λW)−1Xβ+ϵ{\bf y} = ({\bf I} – \lambda {\bf W})^{-1} {\bf X}\beta + \epsilon y=(I–λW)−1Xβ+ϵ

给定X{\bf X}X的值下y{\bf y}y的平均值，被称为y{\bf y}y在X{\bf X}X的条件下的期望。因为 ϵ\epsilonϵ独立于X{\bf X}X， y{\bf y}y在X{\bf X}X下的条件期望是

E(y∣X)=(I–λW)−1XβE({\bf y}|{\bf X}) = ({\bf I} – \lambda {\bf W})^{-1} {\bf X}\beta E(y∣X)=(I–λW)−1Xβ

注意，由于y{\bf y}y是一个向量，这个条件期望公式指定了德克萨斯州每个县的平均值。

我们使用这个等式来定义，当 X{\bf X}X 从一组值变化到另一组值后，y{\bf y}y 的变化效果。具体而言，X0{\bf X_0}X0 为来自原始观察数据中的协变量值；将 X0{\bf X_0}X0 中达拉斯的失业率修改为 10％，其余数据均保持不变，即为 X1{\bf X_1}X1 。通过这种表示，我们计算从 X0{\bf X_0}X0 到 X1{\bf X_1}X1 的改变所引起的德克萨斯州每个县的平均谋杀率的变化。

E(y∣X1)−E(y∣X0)=(I−λW)−1X1β−(I−λW)−1X0β=(I−λW)−1ΔXβ(1)\begin{array}{l}{E\left(\mathbf{y} | \mathbf{X}_{1}\right)-E\left(\mathbf{y} | \mathbf{X}_{0}\right)} \\ {=(\mathbf{I}-\lambda \mathbf{W})^{-1} \mathbf{X}_{1} \beta-(\mathbf{I}-\lambda \mathbf{W})^{-1} \mathbf{X}_{0} \beta} \\ {=(\mathbf{I}-\lambda \mathbf{W})^{-1} \Delta \mathbf{X} \beta}\end{array} \quad{(1)} E(y∣X1)−E(y∣X0)=(I−λW)−1X1β−(I−λW)−1X0β=(I−λW)−1ΔXβ(1)

其中，ΔX=X1–X0\Delta {\bf X}= {\bf X_1} – {\bf X_0}ΔX=X1–X0。

接下来，我们以 SAR 模型为基础，构造用于产生动图的表达式。 SAR 模型之所以被广泛使用，正是因为它们满足稳定条件。这种稳定条件表明逆矩阵 (I–λW)−1({\bf I} – \lambda {\bf W})^{-1}(I–λW)−1 可以写成指数大小逐渐减小的各项之和。这个条件即为：

(I–λW)−1=(I+λW+λ2W2+λ3W3+…)(2)({\bf I} – \lambda {\bf W})^{-1} = ({\bf I} + \lambda {\bf W} + \lambda^2 {\bf W}^2 + \lambda^3 {\bf W}^3 + \ldots) \quad{(2)} (I–λW)−1=(I+λW+λ2W2+λ3W3+…)(2)

将公式 (2) 代入公式 (1)，得到

E(y∣X1)−E(y∣X0)=(I−λW)−1ΔXβ=(I+λW+λ2W2+λ3W3+…)ΔXβ=ΔXβ+λWΔXβ+λ2W2ΔXβ+λ3W3ΔXβ+…(3)\begin{array}{l}{E\left(\mathbf{y} | \mathbf{X}_{1}\right)-E\left(\mathbf{y} | \mathbf{X}_{0}\right)} \\ {=(\mathbf{I}-\lambda \mathbf{W})^{-1} \Delta \mathbf{X} \beta} \\ {=\left(\mathbf{I}+\lambda \mathbf{W}+\lambda^{2} \mathbf{W}^{2}+\lambda^{3} \mathbf{W}^{3}+\ldots\right) \Delta \mathbf{X} \beta} \\ {=\Delta \mathbf{X} \beta+\lambda \mathbf{W} \Delta \mathbf{X} \beta+\lambda^{2} \mathbf{W}^{2} \Delta \mathbf{X} \beta+\lambda^{3} \mathbf{W}^{3} \Delta \mathbf{X} \beta+\ldots}\end{array} \ (3) E(y∣X1)−E(y∣X0)=(I−λW)−1ΔXβ=(I+λW+λ2W2+λ3W3+…)ΔXβ=ΔXβ+λWΔXβ+λ2W2ΔXβ+λ3W3ΔXβ+… (3)

这即为生成动态图形效果的表达式。

公式 (3) 中的每一项都很直观，这在我们所展示的案例中很容易解释：

第一项 (ΔXβ\Delta {\bf X}\betaΔXβ) 是初始效应变化，它只影响达拉斯地区的凶杀率。
第二项 (λWΔXβ\lambda {\bf W} \Delta {\bf X}\betaλWΔXβ) 是达拉斯地区凶杀率的变化对其邻居的影响。
第三项 (λ2W2ΔXβ\lambda^2 {\bf W}^2 \Delta {\bf X}\betaλ2W2ΔXβ) 是达拉斯地区凶杀率的变化对达拉斯邻居的邻居的影响。依次类推，可以得到其他变量的含义。

为溢出效应创建动态图形

我现在描述如何生成动态图。每个图使用公式 (3) 中每一项子集绘制变化图。第一个图仅显示第一个项计算的变化。第二个图仅显示从第一项到第二项计算的变化。第三个图仅显示从第一项到第三项计算的变化，以此类推。

代码的前四个步骤执行以下操作。

计算并绘制 ΔXβ\Delta {\bf X}\betaΔXβ。
计算并绘制 ΔXβ+λWΔXβ\Delta {\bf X} \beta + \lambda {\bf W} \Delta {\bf X}\betaΔXβ+λWΔXβ。
计算并绘制 ΔXβ+λWΔXβ+λ2W2ΔXβ\Delta {\bf X} \beta + \lambda {\bf W} \Delta {\bf X}\beta + \lambda^2 {\bf W}^2 \Delta {\bf X}\betaΔXβ+λWΔXβ+λ2W2ΔXβ。
计算并绘制 ΔXβ+λWΔXβ+λ2W2ΔXβ+λ3W3ΔXβ\Delta {\bf X} \beta + \lambda {\bf W} \Delta {\bf X}\beta + \lambda^2 {\bf W}^2 \Delta {\bf X}\beta + \lambda^3 {\bf W}^3 \Delta {\bf X} \betaΔXβ+λWΔXβ+λ2W2ΔXβ+λ3W3ΔXβ。

步骤 5 到 20 执行类似的操作。最后，结合步骤 1 到步骤 20 中的图形，创建一个动态图形。
下面的代码展示了此过程：

  1 /* get estimate of spatial lag parameter lambda */2 local lambda = _b[W:hrate]34 /* xb based on original data */5 predict xb0, xb67 /* xb based on modified data */8 replace unemployment = 10 if cname == "Dallas"9 predict xb1, xb1011 /* compute the outcome change in the first step */12 generate dy = xb1 - xb013 format dy %9.2f1415 /* Initialize Wy, lamWy, */16 generate Wy = dy17 generate lamWy = dy1819 /* map the outcome change in step 1 */20 grmap dy21 graph export dy_0.png, replace22 local input dy_0.png2324 /* compute the outcome change from step 2 to 11 */25 forvalues p=1/20 {26         spgenerate tmp = W*Wy27         replace lamWy = `lambda'^`p'*tmp28         replace Wy = tmp29         replace dy = dy + lamWy30         grmap dy31         graph export dy_`p'.png, replace32         local input `input' dy_`p'.png33         drop tmp34 }3536 /* convert graphs into a animated graph */37 shell convert -delay 150 -loop 0 `input' glsp.gif3839 /* delete the generated pgn file */40 shell rm -fR *.png

此代码使用由 spregress 命令估计过程中所产生的 ereturn 类的返回值及其相应的 predict命令。

第 2 行将 λ\lambdaλ 的估计值放在局部暂元 lambda 中。
第 5, 7, 8 和 9 分别计算 X0{\bf X_0}X0和X1{\bf X_1}X1 下的 Xβ{\bf X}\betaXβ ，并将其存储在 xb0 和 xb1 中。
第 12 行计算 (ΔXβ\Delta {\bf X}\betaΔXβ) 并将其存储在 dy 中。
第 16 行和第 17 行存储当 p=0p = 0p=0 时 Wpy{\bf W}^{p} {\bf y}Wpy 和 λpWpy\lambda^{p} {\bf W}^{p} {\bf y}λpWpy 的初始值。
第 20-22 行生成动态图中的第一个图。当全部代码完成时，局部暂元 input 将包含用于创建动态图形中的每一张图形。
第 25-34 行计算并绘制剩余项目的图形。第 26 行使用 spgenerate 命令计算 Wpy{\bf W}^{p} {\bf y}Wpy。
第 27-33 行执行与 dy 类似的操作。
在第 37 行，我使用 Linux 工具 “convert” 来组合图形以生成动态图形。在Windows上，我可以使用 FFmpeg 和 Camtasia 等软件。有关详细信息，请参见Chuck Huber的How to create animated graphics using Stata
第 40 行删除所有不必要的 .png 文件。

下面是由此代码创建的动态图形。

小结

在这篇文章中，利用德克萨斯州案例，我讨论SAR模型如何解释溢出效应。我还介绍了如何将溢出效应计算为累加和，并利用累积的总和创建了一个动态图形，说明了这些溢出效应是如何在德克萨斯州的各个县蔓延的。

关于我们

【Stata 连享会(公众号：StataChina)】由中山大学连玉君老师团队创办，旨在定期与大家分享 Stata 应用的各种经验和技巧。
公众号推文同步发布于 CSDN-Stata连享会、简书-Stata连享会和知乎-连玉君Stata专栏。可以在上述网站中搜索关键词Stata或Stata连享会后关注我们。
点击推文底部【阅读原文】可以查看推文中的链接并下载相关资料。
Stata连享会精彩推文1 || 精彩推文2

联系我们

欢迎赐稿： 欢迎将您的文章或笔记投稿至Stata连享会(公众号: StataChina)，我们会保留您的署名；录用稿件达五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。
意见和资料： 欢迎您的宝贵意见，您也可以来信索取推文中提及的程序和数据。
招募英才： 欢迎加入我们的团队，一起学习 Stata。合作编辑或撰写稿件五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。
联系邮件： StataChina@163.com

往期精彩推文

Stata连享会推文列表

Stata连享会精彩推文1 || 精彩推文2