Matlab：莫兰指数的编程实现
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
连享会 · 2022 面板数据因果推断专题
作者：陈子厚 (中共广东省委党校)
邮箱：econometricalc@outlook.com
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
目录
1. 莫兰指数的定义与公式
2. 莫兰指数的 Matlab 编程
2.1 Matlab 的编程实现
2.2 与 Stata 结果的验证
3. 附录
4. 相关推文
计算莫兰指数的软件非常多，例如 Stata，详见连享会推文「Stata：面板数据的莫兰指数计算与散点图绘制-xtmoran」。本文主要借助 Matlab 实现莫兰指数的计算，有编程兴趣的同学，可以将此代码作为练习。
1. 莫兰指数的定义与公式
莫兰指数是全域莫兰指数 (Global Moran's I) 的简称，是澳大利亚统计学家帕特里克·阿尔弗雷德·皮尔斯·莫兰在 1950 年提出的，是最常用的空间自相关指标。其计算公式如下：
公式中 和 分别表示第 个空间单元和第 个空间单元的特征值 (如：人均 GDP)， 为地区个数， 为空间权重矩阵。莫兰指数值在正负 1 之间，正值代表正相关性，负值代表负相关性，值越大表示空间集聚性越强，趋于 0 表示在空间中随机分布。
对于莫兰指数的检验，原假设 为：不存在空间自相关。按照假设检验方法，当区域个数 足够大时，全局莫兰指数近似服从正态分布，即：
莫兰指数统计量的期望公式如下：
莫兰指数统计量的方差公式如下：
2. 莫兰指数的 Matlab 编程
为了方便理解和实现 Matlab 编程，需要将莫兰指数公式简化，并分成两部分：一是莫兰指数值的公式简化；二是莫兰指数的检验公式简化。首先，将莫兰指数的公式简化：
其中，，如果矩阵进行标准化，则 。
其次，将莫兰指数的检验公式简化，其关键在于 的简化，难点在于统计量方差公式中 的简化。这里参考姜磊老师书籍《应用空间计量经济学》(P33 页) 与「Global Moran's I 其他数学计算」，计算如下：
2.1 Matlab 的编程实现
首先，借助上述公式的简化，莫兰指数值的编程实现如下所示。需要注意的是，数据要按列存放，每一列是一年的数据，这里以一列数据 x (即一年数据) 为例。
s0=n;
s=var(x,1);
m=mean(x);
y=0;
for i=1:n
for j=1:n
y=y+w(i,j)*(x(i)-m)*(x(j)-m);
end
end
moran=y/(s*n);
其次，莫兰指数检验的实现。
%% 1. 求解 S1, S2, D
%% 求解 S1
s1=0;
for i=1:n
for j=1:n
s1=s1+(w(i,j)+w(j,i))^2;
end
end
s1=s1/2;
%% 求解 S2
s2=0;
for i=1:n
w12(i)=0;
w21(i)=0;
for j=1:n
w12(i)=w12(i)+w(i,j);
w21(i)=w21(i)+w(j,i);
end
s2=s2+(w12(i)+w21(i))^2;%S2
end
%% 求解 D
d1=x-m;
d2=sum(d1.^4);
d3=(sum(d1.^2))^2;
d=d2/d3;
%% 2. 求解 A, B, C；
A=n*((n^2-3*n+3)*s1-n*s2+3*s0^2);
B=d*n*((n^2-n)*s1-2*n*s2+6*s0^2);
C=s0^2*(n-1)*(n-2)*(n-3);
%% 3. 求解 ei 值, sd 值, Z 值, P 值
ei=-1/(n-1);
ei2=ei^2;
varc=((A-B)/C)-ei2;
sd=varc^(1/2);
z=(moran-ei)/sd;
P = 2 * (1-normcdf(abs(z)));
最后，将结果汇总。
reults=[moran ei sd z P];
以上程序为计算一年的莫兰指数，可以运用循环计算出多年的莫兰指数，演示数据计算的结果如下所示：
莫兰检验结果 (第一列为I(莫兰指数), 第二列为 E(I), 第三列为 SD(I), 第四列为 Z, 第五列为 P)
I =
0.3178 -0.0345 0.1180 2.9847 0.0028
0.3524 -0.0345 0.1148 3.3692 0.0008
0.3526 -0.0345 0.1149 3.3681 0.0008
0.3517 -0.0345 0.1150 3.3583 0.0008
0.3409 -0.0345 0.1154 3.2529 0.0011
0.3475 -0.0345 0.1151 3.3198 0.0009
0.3581 -0.0345 0.1143 3.4358 0.0006
2.2 与 Stata 结果的验证
将数据另存为 dta 格式，并导入 Stata，运用 spatgsa 命令计算莫兰指数。
* 调入空间权重
spatwmat using 矩阵.dta,name(W) standardize
* 调用数据
use 数据.dta,clear
* 全局莫兰指数
spatgsa a2012 a2013 a2014 a2015 a2016 a2017 a2018 ,weights(W) twotail moran
可以看出，Stata 与 Matlab计算的结果是一致的。
Measures of global spatial autocorrelation
Weights matrix
--------------------------------------------------------------
Name: W
Type: Imported (binary)
Row-standardized: Yes
--------------------------------------------------------------
Moran's I
--------------------------------------------------------------
Variables | I E(I) sd(I) z p-value*
--------------------+-----------------------------------------
a2012 | 0.3178 -0.0345 0.1180 2.9847 0.0028
a2013 | 0.3524 -0.0345 0.1148 3.3692 0.0008
a2014 | 0.3526 -0.0345 0.1149 3.3681 0.0008
a2015 | 0.3517 -0.0345 0.1150 3.3583 0.0008
a2016 | 0.3409 -0.0345 0.1154 3.2529 0.0011
a2017 | 0.3475 -0.0345 0.1151 3.3198 0.0009
a2018 | 0.3581 -0.0345 0.1143 3.4358 0.0006
--------------------------------------------------------------
*2-tail test
3. 附录
clear;clc;
A = xlsread('莫兰指数数据.xlsx'); %按列存放, 每一列是一年的数据
W = xlsread('空间邻接矩阵.xlsx');
Y=A(:,3:9); % 按列存放, 每一列是一年的数据
[n,T]=size(Y);
w = normw(W); % 行标准化, 借助 jplv7 工具包的 normw 函数
for k=1:T
x=Y(:,k); % 第 K 年的数据
%% 求莫兰指数值
s0=n;
s=var(x,1);
m=mean(x);
y=0;
for i=1:n
for j=1:n
y=y+w(i,j)*(x(i)-m)*(x(j)-m);
end
end
moran=y/(s*n);
%% 求解 S1
s1=0;
for i=1:n
for j=1:n
s1=s1+(w(i,j)+w(j,i))^2;
end
end
s1=s1/2;
%% 求解 S2
s2=0;
for i=1:n
w12(i)=0;
w21(i)=0;
for j=1:n
w12(i)=w12(i)+w(i,j);
w21(i)=w21(i)+w(j,i);
end
s2=s2+(w12(i)+w21(i))^2;%S2
end
%% 求解 D
d1=x-m;
d2=sum(d1.^4);
d3=(sum(d1.^2))^2;
d=d2/d3;
%% 求解A,B,C
A=n*((n^2-3*n+3)*s1-n*s2+3*s0^2);
B=d*n*((n^2-n)*s1-2*n*s2+6*s0^2);
C=s0^2*(n-1)*(n-2)*(n-3);
%% 求解 ei 值, sd 值, Z 值, P 值
ei=-1/(n-1);
ei2=ei^2;
varc=((A-B)/C)-ei2;
sd=varc^(1/2);
z=(moran-ei)/sd;
P = 2 * (1-normcdf(abs(z)));
reults=[moran ei sd z P];
I(k,:)=reults;
end
% 列表显示莫兰指数结果
fprintf(1,'莫兰检验结果(第一列为I 第二列为E(I) 第三列为SD(I) 第四列为Z 第五列为P)')
I
4. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh 空间计量, m
安装最新版lianxh命令：
ssc install lianxh, replace
专题：Stata命令 Stata：空间计量之用-spmap-绘制地图.md 专题：Stata绘图 Stata空间计量：莫兰指数绘图moranplot命令介绍 专题：空间计量-网络分析 Stata：空间计量模型双权重-spm Stata：批量获取经纬度数据-空间计量 Stata空间计量：STAR-时空自回归模型 Stata：一文遍览Stata官方空间计量命令：sp系列命令 刘迪：Stata空间溢出效应的动态图形-空间计量 空间计量溢出效应的动态GIF演示 空间计量：地理加权归回模型-(GWR)-参数估计 专题：公开课 连享会公开课：空间计量建模与高质量论文撰写--范巧教授 专题：最新课程 ⏩ 2022空间计量专题-连享会 专题：公开课 ⏩重磅推荐：范巧——空间计量经济学公开课 重磅推荐：杨海生——空间计量公开课
课程推荐：面板数据因果推断
主讲老师：徐轶青 (斯坦福大学)
🍓 课程主页：https://gitee.com/arlionn/Course
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。