[转载] 刘迪: Stata空间溢出效应的动态图形 (空间计量)

 

原文链接：/news/0b822f835c9d5.html

原文：How to create animated graphics to illustrate spatial spillover effects中文版PPT 在线浏览6 March 作者：Di Liu, Senior Econometrician

 

 

Stata连享会   主页 || 视频 || 推文

 

连享会 - Stata 暑期班

线上直播 9 天：.7.28-8.7

主讲嘉宾：连玉君 (中山大学) | 江艇 (中国人民大学)

课程主页：/arlionn/PX | 微信版

 

This post shows how to create animated graphics that illustrate the spatial spillover effects generated by a spatial autoregressive (SAR) model. After reading this post, you could create an animated graph like the following.

This post is organized as follows. First, I estimate the parameters of a SAR model. Second, I show why a SAR model can produce spatial spillover effects. Finally, I show how to create an animated graph that illustrates the spatial spillover effects.

A SAR model

I want to analyze the homicide rate in Texas counties as a function of unemployment. I suspect that the homicide rate in one county affects the homicide rate in neighboring counties.

I want to answer two questions.

How can I set up a model that explicitly allows the homicide rate in one county to depend on the homicide rate in neighboring counties?Given my model, if the unemployment rate in Dallas increases to 10%, how would the homicide rate change in the neighboring counties of Dallas ?

Fit a SAR model

A standard linear model for the homicide rate in county iii (hratei)({\bf{hrate}}_{i})(hratei​) as a function of the unemployment rate in that county’s unemploymenti{\bf unemployment}_{i}unemploymenti​ is

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

A SAR model allows (hratei)({\bf hrate}_i)(hratei​) to depend on the homicide rate in neighboring counties. I need some new notation to write down a SAR model. I let (Wi,j)(W_{i,j})(Wi,j​) be a positive number if county jjj is a neighbor of county iii, zero if the jjj is not a neighbor of iii, and zero if j=ij=ij=i, because no county can border itself.

Given this notation, a SAR model that allows the homicide rate in county iii to depend on the homicide rate in neighboring counties can be written as

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​=γ1​j=1∑N​Wi,j​hratej​+β1​unemploymenti​+β0​+ϵi​

where (Wi,j)(W_{i,j})(Wi,j​) defines the closeness between county iii and county jjj. The term ∑j=1NWi,jhratej\sum_{j=1}^N W_{i,j} {\bf hrate}_{j}∑j=1N​Wi,j​hratej​ is a weighted sum of the homicide rates in county iii’s neighboring counties, and it specifies how the homicide rates in neighboring counties affect the homicide rate in county iii.

Stacking the neighborhood information in (Wi,j)(W_{i,j})(Wi,j​) for each county iii produces a matrix W{\bf W}W that records the neighbor information for each county iii. The matrix W{\bf W}W is known as a spatial-weighting matrix.

The spatial-weighting matrix that we are using has a special structure; each element is either a value ccc or zero, where ccc is greater than zero. This type of spatial-weighting matrix is known as a normalized contiguity matrix.

In Stata, we usespmatrixto create a spatial-weighting matrix, and we usespregressto fit a cross-sectional SAR model.

I begin by downloading some data on the homicide rates of U.S. counties from the Stata website and creating a subsample that uses only data on counties in Texas.

. /* Get data for Texas counties' homicide rate */. copy http://www.stata-/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, replacefile texas.dta saved

Intuitively, a file that specifies the borders of all the places of interest is known as a shape filetexas.dtais linked to the Stata version of a shape file that specifies the borders of all the counties in Texas. I now download that dataset from the Stata website and usespsetto show that they are linked.

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

I now usespmatrixto create a normalized contiguity spatial-weighting matrix.

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

Now that I have my data and my spatial-weighting matrix, I can estimate the model parameters.

. /* Estimate SAR model parameters */. spregress hrate unemployment, dvarlag(W) gs2sls(254 observations)(254 observations (places) used)(weighting matrix defines 254 places)Spatial autoregressive modelNumber of obs= 254GS2SLS 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 .1525033.01 0.003.1595237 .7573245_cons | 2.720913 1.6531051.65 0.100 -.5191143 5.960939-------------+----------------------------------------------------------------W |hrate | .3414964 .19148651.78 0.075 -.0338103 .7168031------------------------------------------------------------------------------Wald test of spatial terms:chi2(1) = 3.18 Prob > chi2 = 0.0745

Spatial spillover

Now we are ready to answer the second question. Based on our estimation results fromspregress, we can proceed in three steps.

Predict the homicide rate using original data.Change Dallas’s unemployment rate to 10% and predict the homicide rate pute the difference between two predictions and map it.

. 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 */

The above graph shows that a change in the unemployment rate in Dallas changes the homicide rates in the counties that are near to Dallas, in addition to the homicide rate in Dallas. The change in Dallas spills over to the nearby counties, and the effect is known as a spillover effect.

SAR model and spatial spillover

In this section, I show why a SAR model generates a spillover effect. In the process, I provide a formula for this effect that I use to create the animated graph.

The matrix form for a SAR model is

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

Solving for y{\bf y}y yields

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

The mean value of y{\bf y}y given a value of X{\bf X}X is known as the the expectation of y{\bf y}y conditional on X{\bf X}X. Because ϵ\epsilonϵ is independent of X{\bf X}X , the expectation of y{\bf y}y conditional on X{\bf X}X is

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β

Note that this conditional expectation specifies the mean for each county in Texas because y{\bf y}y is a vector.

We use this equation to define the effect of going from one set of values for X{\bf X}X to another set. In the case at hand, I let X0{\bf X_0}X0​ contain the covariate values in the observed data and let X1{\bf X_1}X1​ contain the same values except that the unemployment rate in Dallas has been set to 10%. With this notation, I see that going from X0{\bf X_0}X0​ to X1{\bf X_1}X1​ causes the mean homicide rates for each county in Texas to change by

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)

where ΔX=X1–X0\Delta {\bf X}= {\bf X_1} – {\bf X_0}ΔX=X1​–X0​.

I now show that a technical condition assumed in SAR models produces an expression for the animated graph. SAR models are widely used because they satisfy a stability condition. Intuitively, this stability condition says that the inverse matrix (I–λW)−1({\bf I} – \lambda {\bf W})^{-1}(I–λW)−1 can be written as a sum of terms that decrease in size exponentially fast. This condition is that

(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)

Plugging the formula from (2) into the effect in (1) yields

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)

which is the expression for the effect that I use to generate the animated graph.

Each term in (3) has some intuition, which is most easily presented in terms of my example. The first term (ΔXβ\Delta {\bf X}\betaΔXβ) is the initial effect of the change, and it affects only the homicide rate in Dallas. The second term (λWΔXβ\lambda {\bf W} \Delta {\bf X}\betaλWΔXβ) is the effect of the change on the outcome in those places that are neighbors of Dallas. The third term (λ2W2ΔXβ\lambda^2 {\bf W}^2 \Delta {\bf X}\betaλ2W2ΔXβ) is the effect of the change on the outcome in those places that are neighbors of neighbors of Dallas. The intuition continues in the pattern for the remaining terms.

Create animated graphs for spillover effects

I now describe how I generate the animated graph. Each graph plots the change using a subset of the terms in (3). The first graph plots the change computed from the first term only. The second graph plots the change computed from the first and second terms only. The third graph plots the change computed from the first three terms only. And so on.

The first four steps of the code do the following.

It computes and plots ΔXβ\Delta {\bf X}\betaΔXβ.It computes and plots ΔXβ+λWΔXβ\Delta {\bf X} \beta + \lambda {\bf W} \Delta {\bf X}\betaΔXβ+λWΔXβ.It compute and plots Δ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β.It computes and plots Δ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β.

Steps 5 through 20 perform the analogous operations.

Finally, combine graphs from step 1 to step 20, and create an animated graph.

Here is the code that implements this process.

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

This code uses the ereturn results produced byspregressabove and its correspondingpredictcommand.

Line 2 puts the estimate of λ\lambdaλ in the local macrolambda.Lines 5, 7, 8, and 9 compute Xβ{\bf X}\betaXβ for X0{\bf X_0}X0​ and X1{\bf X_1}X1​ and store them inxb0andxb1, respectively.Line 12 computes the first term (ΔXβ\Delta {\bf X}\betaΔXβ) and stores it indy.Lines 16 and 17 store the initial values for Wpy{\bf W}^{p} {\bf y}Wpy and λpWpy\lambda^{p} {\bf W}^{p} {\bf y}λpWpy, when p=0p=0p=0.Lines 20–22 produce the first plot in the animated graph. The local macroinputwill contain all the plots used to create the animated graph when the code finishes.Lines 25–34 compute the terms and create the plots for the remaining terms. Line 26 usesspgenerateto compute Wpy{\bf W}^{p} {\bf y}Wpy. Line 27–33 perform operations analogous to those ofdy.In Line 37, I use a Linux tool “convert” to combine the graphs to produce an animated graph. On Windows, I can use software such as FFmpeg and Camtasia. For more details, see How to create animated graphics using Stata by Chuck Huber.Line 40 deletes all the unnecessary.pngfiles.

Here is the animated graph created by this code.

Done and undone

In this post, I discussed spillover effects and why SAR models produce them in the context of an example using the counties in Texas. I also showed how the effects can be computed as an accumulated sum. I used the accumulated sum to create an animated graph that illustrates how the effects spill over in the counties in Texas.

关于我们

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

联系我们

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

往期精彩推文

Stata连享会推文列表

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