Matlab：莫兰指数的编程实现

👇 连享会 · 推文导航 | www.lianxh.cn

• 🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令
• 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享
• 💹 计量：回归分析 | 交乘项-调节  | IV-GMM | 时间序列 | 面板数据  | 空间计量 | Probit-Logit | 分位数回归
• ⛳ 专题：SFA-DEA | 生存分析  | 爬虫 | 机器学习 | 文本分析
• 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching
• 🔨 工具：工具软件 | Markdown | Python-R-Stata
• 🎧 课程：公开课-直播 |  计量专题 | 关于连享会

连享会 · 2022 面板数据因果推断专题

作者：陈子厚 (中共广东省委党校)
邮箱：[e-mail]

温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：

 

---

目录

• 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 和 songbl  GIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
  . ssc install lianxh
  . ssc install songbl
👉  使用：
  . lianxh DID 倍分法
  . songbl all

🍏 关于我们

• 连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。
• 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。