Matlab:莫兰指数的编程实现

发布时间:2022-10-15 阅读 1016

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 系列推文:

PDF下载 - 推文合集

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


目录


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

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

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

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

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

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

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

2. 莫兰指数的 Matlab 编程

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

其中,S0=i=1nj=1nwij,如果矩阵进行标准化,则 S0=n

其次,将莫兰指数的检验公式简化,其关键在于 Var(I) 的简化,难点在于统计量方差公式中 E(I2) 的简化。这里参考姜磊老师书籍《应用空间计量经济学》(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

相关课程

免费公开课

最新课程-直播课

专题 嘉宾 直播/回看视频
最新专题 文本分析、机器学习、效率专题、生存分析等
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 连享会-主页知乎专栏,700+ 推文,实证分析不再抓狂。直播间 有很多视频课程,可以随时观看。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

✏ 连享会-常见问题解答:
https://gitee.com/lianxh/Course/wikis

New! lianxhsongbl 命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh