Stata：如何生成天马行空的空间权重矩阵？

发布时间：2022-09-11 阅读 3945

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下载 - 推文合集

作者：张家星 (中南财经政法大学)
邮箱：zhangjx@stu.zuel.edu.cn

1. 空间权重矩阵简介

空间计量经济学与传统计量经济学的最大区别是引入空间权重矩阵。常用空间权重矩阵包括：邻接矩阵、距离矩阵和嵌套矩阵。在进行空间相关性分析、模型检验以及回归时，选取不同的空间权重矩阵，往往会带来不同的结果。因此，我们有必要去识别不同的空间权重矩阵，掌握其生成原理，以便更好地进行经济学研究。

1.1 邻接矩阵

1.1.1 地理邻接矩阵

W_{i j}^{01} = {\begin{cases} 0 & i 与 j 不相邻 \\ 1 & i 与 j 相邻 \end{cases}

如果地区 $i$ 与地区 $j$ 相邻，则 $W_{i j} = 1$ ；如果地区 $i$ 与地区 $j$ 不相邻，则 $W_{i j} = 0$ 。

1.1.2 广义相邻矩阵

W_{i j}^{s} = {\begin{cases} 0 & i 与 j 不属于同一区域 \\ 1 & i 与 j 属于同一区域 \end{cases}

如果地区 $i$ 与地区 $j$ 属于同一区域，则 $W_{i j} = 1$ ；如果地区 $i$ 与地区 $j$ 不属于同一区域，则 $W_{i j} = 0$ 。

1.2 距离矩阵

1.2.1 地理距离矩阵

反距离矩阵公式：

W_{i j}^{d} = {\begin{cases} \frac{1}{d_{i j}} & i \neq j \\ 0 & i = j \end{cases}

反距离平方矩阵公式：

W_{i j}^{d} = {\begin{cases} \frac{1}{d_{i j}^{2}} & i \neq j \\ 0 & i = j \end{cases}

其中， $d_{i j}$ 为地区 $i$ 和地区 $j$ 之间的距离。

1.2.2 经济距离矩阵

W_{i j}^{e} = {\begin{cases} \frac{1}{| P G D P_{i} - P G D P_{j} |} & i \neq j \\ 0 & i = j \end{cases}

其中， $P G D P_{i}$ 为地区 $i$ 的人均 GDP， $P G D P_{j}$ 为地区 $j$ 的人均 GDP。

1.3 嵌套矩阵

1.3.1 地理经济嵌套矩阵 (相加)

W_{i j}^{d e} = {\begin{cases} α \times W_{i j}^{d} + (1 - α) \times W_{i j}^{e} & i \neq j \\ 0 & i = j \end{cases}

其中， $α$ 为参数， $W_{i j}^{d}$ 为反地理距离矩阵， $W_{i j}^{e}$ 为经济距离矩阵。

1.3.2 地理经济嵌套矩阵 (相乘)

W_{i j}^{d e} = {\begin{cases} W_{i j}^{d} \times W_{i j}^{e} & i \neq j \\ 0 & i = j \end{cases}

其中， $W_{i j}^{d}$ 为反地理距离矩阵， $W_{i j}^{e}$ 为经济距离矩阵。

1.3.3 基于引力模型的矩阵

W_{i j}^{d e} = {\begin{cases} \frac{P G D P_{i} \times P G D P_{j}}{d_{i j}^{2}} & i \neq j \\ 0 & i = j \end{cases}

其中， $P G D P_{i}$ 为地区 $i$ 的人均 GDP， $P G D P_{j}$ 为地区 $j$ 的人均 GDP， $d_{i j}$ 为地区 $i$ 和地区 $j$ 之间的距离。

2. 准备工作

2.1 示例数据

具体示例数据，可通过点击「空间权重矩阵的构建案例数据.zip」下载。主要文件如下：

shapefiles 文件，包括：30省.cpg，30省.dbf，30省.prj，30省.sbn，30省.sbx，30省.shp，30省.shp.xml，30省.shx。8 个文件缺一不可，需全部放在 Stata 当前工作路径之下。
pgdp_15.dta 文件，共包含一个变量。其中 pgdp2015 为 2015 年 30 个省、市、自治区人均 GDP 数据。

2.2 相关命令

在分析过程中，需要用到 spwmatrix、spmatrix、spmat、shp2dta、spshape2dta、mata 等命令，其中 spmatrix 和 spshape2dta 需要在 Stata 15 及以上版本才能运行。相关命令安装如下：

ssc install spwmatrix, replace 
net install st0292.pkg, replace  
ssc install shp2dta, replace

2.3 Mata 语言的基础用法

2.3.1 矩阵常规运算

矩阵乘法运算示意图及代码：

[\begin{matrix} 1 & 2 \\ 3 & 4 \end{matrix}] \times [\begin{matrix} 2 & 3 \\ 4 & 5 \end{matrix}] = [\begin{matrix} 1 \times 2 + 2 \times 4 & 1 \times 3 + 2 \times 5 \\ 3 \times 2 + 4 \times 4 & 3 \times 3 + 4 \times 5 \end{matrix}]

// 矩阵 Wa 和 矩阵 Wb 相乘，得到矩阵 Wab
mata
Wab = Wa * Wb
Wab
end
// mata 和 end 需要同时使用，表示进入 mata 语言

矩阵加法运算示意图及代码：

[\begin{matrix} 1 & 2 \\ 3 & 4 \end{matrix}] + [\begin{matrix} 2 & 3 \\ 4 & 5 \end{matrix}] = [\begin{matrix} 1 + 2 & 2 + 3 \\ 3 + 4 & 4 + 5 \end{matrix}]

// 矩阵 Wa 和 矩阵 Wb 相加，得到矩阵 Wab
mata
Wab = Wa + Wb
Wab
end
// mata 和 end 需要同时使用，表示进入 mata 语言

2.3.2 矩阵元素对元素运算

元素对元素的运算示意图及代码：

[\begin{matrix} 1 & 2 \\ 3 & 4 \end{matrix}] \times [\begin{matrix} 2 & 3 \\ 4 & 5 \end{matrix}] = [\begin{matrix} 1 \times 2 & 2 \times 3 \\ 3 \times 4 & 4 \times 5 \end{matrix}]

两个矩阵各元素对应相乘，即 $W_{i j}^{a} \times W_{i j}^{b}$ 。

// 矩阵 Wa 与矩阵 Wb 中的各元素对应相乘，得到矩阵 Wab
mata
Wab = Wa ：* Wb
Wab
end
// ":" 加上 "+ - * /"，便是元素对元素运算的符号

2.3.3 循环

// for 循环语句
   for (i=1; i<=n; i++) {
      stmts  // stmts 可以理解为循环的具体内容
   }
// 表示 i 从 1 到 n 递增循环

2.3.4 矩阵导入与导出

spmatrix matafromsp 与 spmatrix spfrommata 的基本语法如下：

* 矩阵导入 mata
spmatrix matafromsp matamatrix matavec = spmatname

* 矩阵从 mata 导出
spmatrix spfrommata spmatname = matamatrix matavec, normalize(normalize)

matamatrix：储存在 Mata 中矩阵的名字；
spmatname：储存在 Stata 中矩阵的名字；
matavec：表示唯一识别的一列向量，即观测对象的唯一 id；
normalize(normalize)：标准化选项，包括 none、row、spectral，默认为 spectral 标准化。

2.4 矢量图层获取、编辑与导入

2.4.1 矢量图层获取

标准地图可从「标准地图服务系统网站」下载，不过该网站不提供 SHP 格式的文件，仅能下载 JPG 格式与 EPS 格式的文件，SHP 格式的地图需要通过其他软件转换生成。对于经济学研究者来说，掌握矢量地图的制作方法并不是硬性要求，我们也可通过其他途径获得矢量地图，值得注意的是公开使用地图前需要送到自然资源主管部门审核。

2.4.2 矢量图层编辑

shapefiles 文件的编辑需要使用到 ArcGIS 软件，可网上搜索相关教程。一般来说，只要我们获得所研究区域的 shapefiles 文件，便可在 Stata 中进行如下操作。

2.4.3 矢量图层导入 Stata

将 shapefile 文件导入到 Stata 中的两个常用命令分别是 shp2dta 和 spshape2dta，后者为 Stata 空间计量官方命令。

. shp2dta using 30省.shp, database(province30_dbf) coordinates(province30_shp) ///
>     genid(geoid) gencentroids(zb) replace 

. use province30_dbf.dta, clear
. list 省 x_zb y_zb in 1/5

     +------------------------------------+
     |           省       x_zb       y_zb |
     |------------------------------------|
  1. |       北京市   116.4123   40.18554 |
  2. |       天津市   117.3341   39.29343 |
  3. |       河北省   116.1374   39.54777 |
  4. |       山西省   112.2895   37.57175 |
  5. |  内蒙古自治区   113.9221   44.08926 |
     +------------------------------------+

shp2dta 命令会新生成 province30_dbf.dta 和 province30_shp.dta 两个文件。在上述 shp2dta 命令语句中，各部分含义如下：

30省.shp：30省 shapefiles 文件名；
database(province30_dbf)：新生成的包含 30省.dbf 中信息的 .dta 文件，括号内为自定义名字；
coordinates(province30_shp)：新生成的包含 30省.shp 中信息的 .dta 文件，括号内为自定义名字；
genid(geoid)：为备选选项，生成一个唯一 id，括号内为自定义名字；
gencentroids(zb)：为备选选项，生成质心 x、y 坐标，括号内为自定义名字；

注意：对于距离矩阵或者经济距离矩阵来说，若 .dta 文件中已包含经纬度坐标和唯一 id ，则不需要 shapefiles 文件去生成矩阵。但是对于邻接矩阵来说，则必须需要 shapefiles 文件才能生成矩阵。因为，仅有经纬度数据是没办法判断两个区域是否相邻。

. spshape2dta 30省.shp, replace
. use 30省.dta, clear
. list 省 _CX _CY in 1/5

     +--------------------------------------+
     |           省         _CX         _CY |
     |--------------------------------------|
  1. |       北京市   116.41228   40.185543 |
  2. |       天津市   117.33404   39.293426 |
  3. |       河北省    116.1374    39.54777 |
  4. |       山西省   112.28951   37.571746 |
  5. |  内蒙古自治区    113.9221   44.089265 |
     +--------------------------------------+

spshape2dta 会自动生成所需要的唯一 id、经纬度坐标，同时会自动生成两个新文件 30省.dta 和 30省_shp.dta。在上述 spshape2dta 命令语句中，各部分含义如下：

30省.dta：新生成的包含 30省.dbf 中信息的 .dta 文件，命令自动生成的名字；
30省_shp.dta：新生成的包含 30省.shp 中信息的 .dta 文件，命令自动生成的名字；
CX 和 _CY：分别为经纬度质心坐标。

shp2dta 和 spshape2dta 均可将 shapefiles 文件导入为 Stata 可识别的 .dta 文件，二者并无明显差别。如果打算采用 Stata 空间计量官方系列命令 SP 生成空间邻接矩阵并进行后续的回归分析，建议采用 spshape2dta。后文主要使用 shp2dta 生成的 province30_dbf.dta 和 province30_shp.dta 来产生空间邻接矩阵。

3. 空间权重矩阵的构建

3.1 地理邻接矩阵

具体公式如下：

W_{i j}^{01} = {\begin{cases} 0 & i 与 j 不相邻 \\ 1 & i 与 j 相邻 \end{cases}

3.1.1 使用 spmat 生成地理邻接矩阵

第一步：导入数据生成邻接矩阵

. use province30_dbf.dta, clear
. spmat contiguity W01 using province30_shp,  ///
>     id(geoid) rook normalize(row) replace    
warning: spatial-weighting matrix contains 1 island

spmat 命令语句含义如下：

contiguity：生成邻接矩阵的选项；
W01：新生成的矩阵名字；
province30_shp：2.43 中生成的包含 30 省地理信息的 .dta 文件；
id(geoid)：2.43 中生成的唯一 id；
rook：默认为 queen 邻接，加 rook 选项后为 rook 邻接。其中，rook 邻接是两个区域有一个面邻接则为 1，queen 邻接是两个区域有一个点邻接则为 1；
normalize(row)：行标准化，还有 minmax 和 spectral 两个选项。

注意：新生成的矩阵中存在孤岛，即有一个观测值周围没有邻居，这就会导致矩阵中有一行元素全部为 0，会影响到后续的空间相关性检验与空间回归，所以我们需要手动或自动为它找一个邻居，这也是生成地理邻接矩阵最容易被大家忽视的一步。

第二步：将矩阵保存成 .spmat 文件

. spmat save W01 using W01.spmat, replace

第三步：将矩阵保存成 .dta 文件

常用空间命令当中所使用的矩阵均大多是 .dta 文件，所以本步是将 W01.spmat 转化为 W01.dta。

. clear all
. spmat use W01 using W01.spmat , replace
. spmat getmatrix W01 W
. getmata (x*)=W
. save W01.dta ,replace

注意：x*，即 x1-xn，是按照原来唯一 id 即 (geoid) 顺序排列的。

第四步：检查是否存在孤岛

Stata 中提示：warning: spatial-weighting matrix contains 1 island，表明新生成的矩阵中存在孤岛，所以我们需要手动或自动为它找一个邻居。本文采取手动方式，感兴趣的朋友也可 help spwmatfill 研究如何自动生成。

注意：spwmatfill 仅适合由 GeoDa 软件生成的 .gal 格式的空间邻接矩阵。

. * 海南 id：21 是孤岛，根据文献手动为其添加一个邻居 广东 id：19
. replace x19 = 1 in 21
. * 如果广东成为海南的邻居，那么广东的邻居数也应该由 4 个增加至 5 个
. replace x13 = 0.2 in 19
. replace x14 = 0.2 in 19
. replace x18 = 0.2 in 19
. replace x20 = 0.2 in 19
. replace x21 = 0.2 in 19
. save W01.dta, replace

至此，我们成功生成出采用 rook 邻接并经过行标准化后的地理邻接矩阵 W01。

3.1.2 其他生成地理邻接矩阵的命令

除 spmat 外，spmatrix 与 spwmatrix 也可生成地理邻接矩阵，矩阵保存及导出方法类似，故不再详细说明，本文仅列出这两个命令的核心用法及注意事项。与此同时，本文提供一种导入 GeoDa 软件所生成的权重矩阵的方法。

spmatrix 生成地理邻接矩阵的核心代码：

. use 30省.dta, clear
. spset //告诉 stata 是空间数据，使用 spmatrix 必要步骤
. spset, modify coordsys(latlong, kilometers) //修改单位参考系
. spmatrix create contiguity W01_2,normalize(row) replace
. spmatrix export W01_2 using W01_2.txt       //导出矩阵

其中，W01_2 新生成的自定义矩阵名字；normalize(row) 中 row 为行标准化，默认为 spectral 标准化，其他标准化选项 minmax 和 none。

spwmatrix 生成地理邻接矩阵的核心代码：

spwmatrix gecon varlist [if] [in], wname(wght_name) ///
   wtype(bin) cart r(#) dband(numlist)              ///
   alpha(#) knn(#) econvar(varname1) beta(#) Other_options

其中，wtype(bin) 中的 bin 代表生成邻接矩阵的选项。

spwmatrix 导入 GeoDa 软件所生成的地理邻接矩阵的核心代码：

spwmatrix import using 30省.gal,wname(W01_3) row ///
   xport(W01_3, dat) replace
svmat W01_3
matname W01_3 x*
save W01_3.dta

其中，

30省.gal：GeoDa 软件所生成的地理邻接矩阵文件；
wname(W01_3)：W01_3 为自定义名字；
row：行标准化；
svmat：将矩阵导入到变量中，详见 help 文件；
matname：为矩阵的行列改名。

3.2 广义相邻矩阵

广义相邻矩阵的一个例子：全国30个省、市、自治区可划分东中西三大区域，生成一列地区虚拟变量 region，其中东部为 1，中部为 2，西部为 3。如果两个省属于相同区域 (即 region 相同) 则 $W_{i j}$ 为1，如果两个省份分属不同区域 (即 region 不同) 则 $W_{i j}$ 为 0。具体公式如下：

W_{i j}^{s} = {\begin{cases} 0 & i 与 j 不属于同一区域 \\ 1 & i 与 j 属于同一区域 \end{cases}

. use province30_dbf.dta,clear   
. * 生成东中西虚拟变量，其中 1 为东部，2 为中部，3 为西部   
. gen region = 1
. replace region=2 if geoid==4 | geoid==7 | geoid==8  | geoid==12|    ///
>                     geoid==14 | geoid==16 | geoid==17 | geoid==18       
. replace region=3 if geoid==5 | geoid==20 | geoid==22  | geoid==23|  ///
>                     geoid==24 | geoid==25 | geoid==26 | geoid==27|  ///
>                     geoid==28 | geoid==29 | geoid==30
. save province30_dbf_region.dta,replace

. * 生成广义相邻矩阵
. use province30_dbf_region.dta, clear
. spwmatrix socio region, wname(Wsoc01) wtype(socnet)  ///
>           idvar(geoid) xport(Wsoc01,txt) replace

. clear all
. spmat import Wsoc01 using "Wsoc01.txt", replace
. spmat getmatrix Wsoc01 W
. getmata (x*)=W
. save Wsoc01.dta ,replace

其中，

spwmatrix socio：社会网络矩阵生成命令；
region：province30_dbf_region.dta 中一个虚拟变量，东部为 1，中部为 2，西部为 3；
wname(Wsoc01)：自定义的矩阵名字；
wtype(socnet)：生成矩阵的类型选社会网络矩阵；
idvar(id)：唯一 id；
xport(Wsoc01,txt)：导出成 .txt 文件，.txt 文件最大的好处就在于几乎每一个空间计量命令都可以识别，从而可以实现跨命令操作；
getmata：可将矩阵导入到数据中，从而生成相应的 .dta 文件。

3.3 地理距离矩阵

Stata 生成地理距离矩阵的命令主要有 spwmatrix 与 spmatrix，同时也可以采用 spmat 命令导入 GeoDa 软件所生成的地理距离矩阵。本文主要基于 spwmatrix 生成地理距离矩阵、地理距离平方矩阵，并保存成相应格式文件。

3.3.1 spwmatrix 生成地理距离矩阵

对于反距离矩阵与反距离平方矩阵来说，最容易被忽视的就是距离的单位。kilometers 和 miles 会产生不同的结果，spwmatrix 默认的距离单位是 kilometers。

反地理距离矩阵具体公式如下：

W_{i j}^{d} = {\begin{cases} \frac{1}{d_{i j}} & i \neq j \\ 0 & i = j \end{cases}

. clear all
. use province30_dbf.dta,clear
. spwmatrix gecon  y_zb x_zb , wn(Wd) wtype(inv) ///
>     alpha(1) xport(Wd,txt) replace             // 加上 row，即可对矩阵行标准化
. clear all
. spmat import Wd using "Wd.txt", replace
. spmat getmatrix Wd W
. getmata (x*)=W
. save Wd.dta, replace
. use Wd.dta, clear

其中，

province30_dbf.dta：2.43 中生成的包含原 shapefiles 文件图层信息的 .dta 文件
y_zb x_zb：经度、纬度，一定要按照这个顺序排列；
wn(Wd)：新生成的矩阵的名字；
wtype(inv)：inv 选项为距离衰减矩阵 (反距离矩阵)；
alpha(1)：默认为 1，若为 2 则生成的是反距离平方矩阵；
xport(Wd,txt)：生成矩阵的同时，导出成 .txt 文件，生成嵌套矩阵时会有大用处；
row：加上则为行标准化。

反地理距离平方矩阵具体公式如下：

W_{i j}^{d} = {\begin{cases} \frac{1}{d_{i j}^{2}} & i \neq j \\ 0 & i = j \end{cases}

. use province30_dbf.dta, clear
. spwmatrix gecon  y_zb x_zb , wn(Wd2) wtype(inv) ///
>     alpha(2) xport(Wd2,txt) replace              // 加上 row，即可对矩阵行标准化
. clear all
. spmat import Wd2 using "Wd2.txt", replace
. spmat getmatrix Wd2 W
. getmata (x*)=W
. save Wd2.dta, replace
. use Wd2, clear

3.3.2 其他生成地理距离矩阵的命令

spmatrix 命令生成反距离矩阵：

. clear all
. spshape2dta 30省.shp, replace 
. use 30省.dta, clear
. spset, modify coordsys(latlong, kilometers)            //更改单位成 kilometers
. spmatrix create idistance Wind,normalize(none) replace //未经过标准化     
. spmatrix export Wind using Wind.txt, replace           //导出成 .txt 文件

. clear all
. spmat import Wind using "Wind.txt", replace
. spmat getmatrix Wind W
. getmata (x*)=W
. save Wind.dta, replace

其中，

spset：将数据定义为空间数据
modify coordsys(latlong, kilometers)：将单位更改为 kilometers，默认单位为 plana，备选中还有 miles；
spmatrix create idistance：生成反距离矩阵的命令；
normalize(none)：是否进行标准化选项，如行标准化，则替换为 row。

注：spmatrix 无法直接生成反距离平方矩阵，spwmatrix 与 spmatrix 两种命令生成的反距离矩阵近似相等 (小数点后三位保持一致)。

导入 GeoDa 软件所生成的反距离距离矩阵：

spmat import A using 30省.gwt, geoda replace //导入 geoda 生成的距离矩阵，命名为 A
spmat getmatrix A W
getmata (x*)=W
save 30省.dta, replace

其中，

30省.gwt：GeoDa 软件生成的反距离矩阵文件；
GeoDa软件：可在软件内自行选取反距离矩阵或反距离平方矩阵，同时可以手动更改距离单位。

注：第二行与第三行代码连用可以生成 x1，x2...x30 存放着刚刚导入进来的反距离矩阵，x1，x2，...，x30 的排列顺序与 shapefiles 文件生成矩阵时所用的唯一 id 顺序一样。

3.4 经济距离矩阵

经济距离的形式有很多种，本文以经济变量人均 GDP 构造经济距离矩阵，具体公式如下：

W_{i j}^{e} = {\begin{cases} \frac{1}{| P G D P_{i} - P G D P_{j} |} & i \neq j \\ 0 & i = j \end{cases}

除此之外，也可以通过更换经济变量或改变计算公式来构造适合您研究的经济距离矩阵，本文以常见的人均 GDP 距离矩阵为例，提供生成经济距离矩阵的代码。在此代码之下，稍作调整，便可生成其他类型的经济距离矩阵。

第一步：导入、合并经济数据

. use province30_dbf.dta ,clear 
. rename 省 province
. merge 1:1 province using "pgdp_15.dta"
. sort geoid  //按照 geoid 排序
. drop _merge //删掉多余的值
. save province30_pgdp15.dta, replace

其中，

province30_dbf.dta：2.43 中生成的带有经纬度坐标的文件；
pgdp_15.dta：带有 2015 年人均 GDP 的文件；
merge：将两个文件合并。

第二步：生成经济距离矩阵

本文是以 2015 年的人均 GDP 作为经济变量，考虑到在后续空间相关性分析及空间回归等命令中所使用的矩阵均为截面矩阵，只取一年经济数据或有不妥，所以文献中常用的一个办法就是对若干年份的人均 GDP 求平均值。本文旨在提供方法，故仅以 2015 年的人均 GDP 为例生成经济距离矩阵。

. clear all
. use province30_pgdp15.dta,clear 
. tomata pgdp2015 
. mata
---------------------- mata (type end to exit) -------------------
:     W = J(30,30,0)
:     id = 1::30
:     for (i=1; i<=length(pgdp2015); i++) {
>         for (j=1; j<=length(pgdp2015); j++) {
>             wi=pgdp2015[i]; wj=pgdp2015[j]
>             if (i!=j)
>             W[i,j]=1/abs(wi-wj)
>         }                
>     }
: end    
------------------------------------------------------------------

其中，

tomata pgdp2015：将 pgdp2015 导入 Mata 中；
W = J(30,30,0)：生成 30 行 30 列，元素全部为 0 的矩阵；
id = 1::30：生成一列 id；
i<=length(pgdp2015);i++)：length() 表示循环次数，本文以变量 pgdp2015 的观测数作为循环次数，故表示为 i 从 1 到 30 递增循环；
if (i!=j)：如果 $i$ 与 $j$ 不相等，便进行下述循环；
W[i,j]=1/abs(wi-wj)：矩阵计算公式的 Stata 表达。

第三步：导出生成的矩阵到 Stata 中并保存成相应文件

. spmatrix spfrommata We_2015 = W id, normalize(none) 
. spmatrix dir

------------------------------------------------------------------------
   Weighting matrix name           N x N      Type         Normalization
------------------------------------------------------------------------
                 We_2015         30 x 30    custom           none
------------------------------------------------------------------------

. spmatrix export We_2015 using We_2015.txt, replace
. clear
. getmata (x*)=W
. save We_2015.dta, replace

其中，

spmatrix spfrommata：从 Mata 中将矩阵导入；
stataspmatname：新生成 Stata 矩阵的自定义名字，如 We_2015；
matamatrix：原始 Mata 矩阵的名字，如 W；
normalize(none)：标准化选项，none 为不标准化，row 为行标准化。

注：这里之所以不对新生成的矩阵进行行标准化，是因为后续地理经济嵌套矩阵的生成会使用到这个矩阵。根据您的研究需要，可自由调整标准化选项。

3.5 地理经济嵌套矩阵 (相加)

经过上述几步，我们已经生成了常见的三种空间权重矩阵，分别是地理邻接矩阵 W01，反地理距离矩阵 Wd，反地理距离平方矩阵 Wd2，经济地理矩阵 We_2015，以及广义相邻矩阵 Wsoc01。

从这一小节开始，我们便以「2.3 Mata 语言的基础用法」中提供的代码为基础，进行矩阵运算，以期生成我们想要的空间权重矩阵。.txt 文件能够在不同空间命令中自由穿梭，故在使用上述命令生成空间权重矩阵过程中，均附带导出生成 .txt 文件。

第一步：地理经济嵌套矩阵 (相加) 生成

W_{i j}^{d e} = {\begin{cases} α \times W_{i j}^{d} + (1 - α) \times W_{i j}^{e} & i \neq j \\ 0 & i = j \end{cases}

为了简便，本文将地理距离矩阵和经济距离矩阵的权重均设置为 0.5，根据研究需要，您可设置不同比例权重。

. spmatrix import Wd using Wd.txt ,replace
. spmatrix import We using We_2015.txt ,replace
. spmatrix matafromsp W1 id = Wd
. spmatrix matafromsp W2 id = We
. mata
---------------------- mata (type end to exit) ---------------------
:     W12 = 0.5:*W1 :+ 0.5:*W2
: end
--------------------------------------------------------------------
. spmatrix spfrommata Wde = W12 id, normalize(none)

思路：首先，借助 spmatrix 导入两个 .txt 格式的矩阵，分别是上几步所生成的反距离矩阵以及经济距离矩阵；其次，通过 spmatrix matafromsp 将矩阵导入 Mata 中；然后，在 Mata 中进行矩阵运算；最后，从 Mata 中，将计算好的矩阵导出到 stata 中。

注：以上矩阵运算均为元素对元素运算 (element_by_element)。

第二步：删除多余缓存矩阵

. mata: mata drop W12 W1 W2 //删掉 Mata 中暂存的矩阵
. spmatrix drop  Wd
. spmatrix drop  We         //删掉 spmatrix 中暂存的矩阵
. spmatrix dir              //矩阵展示   

------------------------------------------------------------------------
   Weighting matrix name           N x N      Type         Normalization
------------------------------------------------------------------------
                     Wde         30 x 30    custom           none
                 We_2015         30 x 30    custom           none
------------------------------------------------------------------------

第三步：导出、保存矩阵

. spmatrix export Wde using Wde_2015.txt,replace
. clear all
. spmat import Wde_2015 using "Wde_2015.txt", replace
. spmat getmatrix Wde_2015 W
. getmata (x*)=W
. save Wde_2015.dta, replace

需要注意一点的是，凡是涉及到经济变量的矩阵，到底是选取某一年的经济变量，还是选取多年经济变量的平均值，这些细节，需要您根据研究内容做出适当选择。不过生成权重矩阵的方法是一致的。

3.6 地理经济嵌套矩阵 (相乘)

具体公式如下：

W_{i j}^{d e} = {\begin{cases} W_{i j}^{d} \times W_{i j}^{e} & i \neq j \\ 0 & i = j \end{cases}

3.6.1 采用 Mata 手动生成地理经济嵌套矩阵 (相乘)

. * 生成矩阵
. spmatrix import Wd using Wd.txt, replace
. spmatrix import We using We_2015.txt, replace
. spmatrix matafromsp W1 id = Wd
. spmatrix matafromsp W2 id = We
. mata
------------------ mata (type end to exit) --------------------
:     W12 = W1 :* W2
: end
---------------------------------------------------------------

. * 删除多余矩阵
. spmatrix spfrommata Wde = W12 id, normalize(none)       
. mata: mata drop W12 W1 W2 //删掉 Mata 中暂存的矩阵
. spmatrix drop  Wd
. spmatrix drop  We         //删掉 spmatrix 中暂存的矩阵
. spmatrix dir              //矩阵展示

------------------------------------------------------------------------
   Weighting matrix name           N x N      Type         Normalization
------------------------------------------------------------------------
                     Wde         30 x 30    custom           none
------------------------------------------------------------------------

. * 导出矩阵
. spmatrix export Wde using Wde_2015_2.txt, replace
. clear all
. spmat import Wde_2015_2 using "Wde_2015_2.txt", replace
. spmat getmatrix Wde_2015_2 W
. getmata (x*)=W
. save Wde_2015_2.dta, replace

3.6.2 采用 spwmatrix 自动生成地理经济嵌套矩阵 (相乘)

spwmatrix 所采用的计算公式如下：

W_{i j}^{d e} = {\begin{cases} \frac{1}{| P G D P_{i} - P G D P_{j} + 1 |} \times e^{- β \times d_{i j}} & i \neq j \\ 0 & i = j \end{cases}

$β$ 默认为 1， $P G D P_{i}$ 为地区 i 的人均 GDP， $P G D P_{j}$ 为地区 $j$ 的人均 GDP， $d_{i j}$ 为地区 $i$ 和地区 $j$ 之间的距离。

. * 生成矩阵
. clear all
. use province30_pgdp15.dta,clear
. spwmatrix gecon  x_zb y_zb, wn(Wecon_2015) wtype(invecon) cart ///
>     econvar(pgdp2015) xport(Wecon_2015, txt) replace

. * 导出矩阵
. clear al
. spmat import Wecon_2015 using "Wecon_2015.txt", replace
. spmat getmatrix Wecon_2015 W
. getmata (x*)=W
. save Wecon_2015.dta, replace
. use Wecon_2015, clear

其中，

province30_pgdp15.dta：带有经纬度及经济变量数据的文件；
spwmatrix gecon-：矩阵生成命令；
x_zb y_zb：经纬度坐标；
wn(Wecon_2015)：自定义矩阵名字；
wtype(invecon)：选择生成反经济距离选项；
cart：表示经纬度坐标采用笛卡尔坐标；
econvar(pgdp2015)：声明 pgdp2015 这一列代表经济变量；
xport(Wecon_2015,txt)：导出成 .txt 文件。

3.7 传统引力模型矩阵

具体公式如下：

W_{i j}^{d e} = {\begin{cases} \frac{P G D P_{i} \times P G D P_{j}}{d_{i j}^{2}} & i \neq j \\ 0 & i = j \end{cases}

$P G D P_{i}$ 为地区 $i$ 的人均 GDP， $P G D P_{j}$ 为地区 $j$ 的人均 GDP， $d_{i j}$ 为地区 $i$ 和地区 $j$ 之间的距离。

思路：分别生成反距离平方矩阵和人均 GDP 两两相乘矩阵，然后两个矩阵再相乘。

生成人均 GDP 两两相乘矩阵：

. * 矩阵计算
. clear all
. use province30_pgdp15.dta,clear 
. tomata pgdp2015 
. mata
---------------------- mata (type end to exit) ------------------
:     W = J(30,30,0)
:     id = 1::30
:     for (i=1; i<=length(pgdp2015);i++) {
>         for (j=1; j<=length(pgdp2015);j++) {
>             Wi=pgdp2015[i]; Wj=pgdp2015[j]
>                 if (i!=j)
>             W[i,j]=Wi*Wj       
>         }                
>     }
: end 
----------------------------------------------------------------

. * 将 Mata 生成的矩阵导入 Stata 中
. spmatrix spfrommata We2_2015 = W id ,normalize(none) 
. spmatrix dir
------------------------------------------------------------------------
   Weighting matrix name           N x N      Type         Normalization
------------------------------------------------------------------------
                We2_2015         30 x 30    custom           none
------------------------------------------------------------------------

. spmatrix export We2_2015 using We2_2015.txt , replace
 
. * 导出成 dta 文件  
. clear
. getmata (x*)=W
. save We2_2015.dta, replace
. use We2_2015.dta, clear

生成传统引力模型矩阵：

. spmatrix import Wd using Wd2.txt ,replace
. spmatrix import We using We_2015.txt ,replace
. spmatrix matafromsp W1 id = Wd
. spmatrix matafromsp W2 id = We

. mata
-------------------------- mata (type end to exit) --------------------
: W12 = W1 :* W2
: end
-----------------------------------------------------------------------

. spmatrix spfrommata Wde = W12 id, normalize(none)
. mata: mata drop W12 W1 W2 //删掉 Mata 中暂存的矩阵
. spmatrix drop  Wd
. spmatrix drop  We         //删掉 spmatrix 中暂存的矩阵
. spmatrix dir              //矩阵展示

------------------------------------------------------------------------
   Weighting matrix name           N x N      Type         Normalization
------------------------------------------------------------------------
                     Wde         30 x 30    custom           none
                We2_2015         30 x 30    custom           none
------------------------------------------------------------------------

. spmatrix export Wde using Wde2_2015.txt,replace
. clear all
. spmat import Wde2_2015 using "Wde2_2015.txt", replace
. spmat getmatrix Wde2_2015 W
. getmata (x*)=W
. save Wde2_2015.dta, replace
. use Wde2_2015, clear

4. 结语

至此，总结一下本文所生成的空间权重矩阵：

地理邻接矩阵：W01
广义”相邻“矩阵：Wsoc01
反地理距离矩阵：Wd
反地理距离平方矩阵：Wd2
经济地理矩阵：We_2015
地理经济嵌套矩阵 (相加)：Wde_2015
地理经济嵌套矩阵 (相乘)：Wde_2015_2 和 Wecon_2015
引力模型矩阵：Wde2_2015

在后续空间相关性分析、模型检验及回归分析中，不同命令对于空间权重矩阵格式的要求不尽相同，不过大多数命令都可以支持识别 .txt 文件。

5. 相关推文

Note：产生如下推文列表的 Stata 命令为：
lianxh 空间计量, m
安装最新版 lianxh 命令：
ssc install lianxh, replace

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

Stata：如何生成天马行空的空间权重矩阵？

1. 空间权重矩阵简介

1.1 邻接矩阵

1.1.1 地理邻接矩阵

1.1.2 广义相邻矩阵

1.2 距离矩阵

1.2.1 地理距离矩阵

1.2.2 经济距离矩阵

1.3 嵌套矩阵

1.3.1 地理经济嵌套矩阵 (相加)

1.3.2 地理经济嵌套矩阵 (相乘)

1.3.3 基于引力模型的矩阵

2. 准备工作

2.1 示例数据

2.2 相关命令

2.3 Mata 语言的基础用法

2.3.1 矩阵常规运算

2.3.2 矩阵元素对元素运算

2.3.3 循环

2.3.4 矩阵导入与导出

2.4 矢量图层获取、编辑与导入

2.4.1 矢量图层获取

2.4.2 矢量图层编辑

2.4.3 矢量图层导入 Stata

3. 空间权重矩阵的构建

3.1 地理邻接矩阵

3.1.1 使用 spmat 生成地理邻接矩阵

3.1.2 其他生成地理邻接矩阵的命令

3.2 广义相邻矩阵

3.3 地理距离矩阵

3.3.1 spwmatrix 生成地理距离矩阵

3.3.2 其他生成地理距离矩阵的命令

3.4 经济距离矩阵

3.5 地理经济嵌套矩阵 (相加)

3.6 地理经济嵌套矩阵 (相乘)

3.6.1 采用 Mata 手动生成地理经济嵌套矩阵 (相乘)

3.6.2 采用 spwmatrix 自动生成地理经济嵌套矩阵 (相乘)

3.7 传统引力模型矩阵

4. 结语

5. 相关推文

相关课程

免费公开课

最新课程-直播课

⛳ 课程主页

⛳ 课程主页

关于我们

热门资讯