Matlab：莫兰指数的编程实现

发布时间：2022-10-15 阅读 1498

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

作者：陈子厚 (中共广东省委党校)
邮箱：econometricalc@outlook.com

计算莫兰指数的软件非常多，例如 Stata，详见连享会推文「Stata：面板数据的莫兰指数计算与散点图绘制-xtmoran」。本文主要借助 Matlab 实现莫兰指数的计算，有编程兴趣的同学，可以将此代码作为练习。

1. 莫兰指数的定义与公式

莫兰指数是全域莫兰指数 (Global Moran's I) 的简称，是澳大利亚统计学家帕特里克·阿尔弗雷德·皮尔斯·莫兰在 1950 年提出的，是最常用的空间自相关指标。其计算公式如下：

I = \frac{n}{\sum_{i = 1}^{n} \sum_{j = 1, j \neq i}^{n} W_{i j}} \cdot \frac{\sum_{i = 1}^{n} \sum_{j = 1}^{n} W_{i j} (x_{i} - \bar{x}) (x_{j} - \bar{x})}{\sum_{i = 1}^{n} {(x_{i} - \bar{x})}^{2}}

公式中 $x_{i}$ 和 $x_{j}$ 分别表示第 $i$ 个空间单元和第 $j$ 个空间单元的特征值 (如：人均 GDP)， $n$ 为地区个数， $W_{i j}$ 为空间权重矩阵。莫兰指数值在正负 1 之间，正值代表正相关性，负值代表负相关性，值越大表示空间集聚性越强，趋于 0 表示在空间中随机分布。

对于莫兰指数的检验，原假设 $H_{0}$ 为：不存在空间自相关。按照假设检验方法，当区域个数 $n$ 足够大时，全局莫兰指数近似服从正态分布，即：

Z = \frac{I - E (I)}{\sqrt{var (I)}}

莫兰指数统计量的期望公式如下：

E (I) = - \frac{1}{n - 1}

莫兰指数统计量的方差公式如下：

Var (I) = E (I^{2}) - E (I)^{2}

2. 莫兰指数的 Matlab 编程

为了方便理解和实现 Matlab 编程，需要将莫兰指数公式简化，并分成两部分：一是莫兰指数值的公式简化；二是莫兰指数的检验公式简化。首先，将莫兰指数的公式简化：

\begin{aligned} I = & \frac{N}{S_{0}} \cdot \frac{\sum_{i}^{N} \sum_{j}^{N} w_{i j} (x_{i} - \bar{x}) (x_{j} - \bar{x})}{\sum_{i}^{n} {(x_{i} - \bar{x})}^{2}} \\ = & \frac{\sum_{i}^{N} \sum_{j}^{N} W_{i j} (x_{i} - \bar{x}) (x_{j} - \bar{x})}{S_{0} \cdot Var (x)} \end{aligned}

其中， $S_{0} = \sum_{i = 1}^{n} \sum_{j = 1}^{n} w_{i j}$ ，如果矩阵进行标准化，则 $S_{0} = n$ 。

其次，将莫兰指数的检验公式简化，其关键在于 $V a r (I)$ 的简化，难点在于统计量方差公式中 $E (I^{2})$ 的简化。这里参考姜磊老师书籍《应用空间计量经济学》(P33 页) 与「Global Moran's I 其他数学计算」，计算如下：

\begin{aligned} E (I^{2}) & = \frac{A - B}{C} \\ A & = n ((n^{2} - 3 n + 3) S_{1} - n S_{2} + 3 S_{0}^{2}) \\ B & = D ((n^{2} - n) S_{1} - 2 n S_{2} + 6 S_{0}^{2}) \\ C & = (n - 1) (n - 2) (n - 3) S_{0}^{2} \\ D & = \frac{\sum_{i = 1}^{n} (x_{i} - \bar{x})^{4}}{{(\sum_{i = 1}^{n} (x_{i} - \bar{x})^{2})}^{2}} \\ S_{1} & = (1 / 2) \sum_{i = 1}^{n} \sum_{j = 1}^{n} {(w_{i, j} + w_{j, i})}^{2} \\ S_{2} & = \sum_{i = 1}^{n} {(\sum_{j = 1}^{n} w_{i, j} + \sum_{j = 1}^{n} w_{j, i})}^{2} \end{aligned}

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

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

Matlab：莫兰指数的编程实现

1. 莫兰指数的定义与公式

2. 莫兰指数的 Matlab 编程

2.1 Matlab 的编程实现

2.2 与 Stata 结果的验证

3. 附录

4. 相关推文

相关课程

免费公开课

最新课程-直播课

⛳ 课程主页

⛳ 课程主页

关于我们

热门资讯