Stata连享会 主页 || 视频 || 推文 || 知乎 || Bilibili 站
温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc
,ihelp
,rdbalance
,gitee
,installpkg
⛳ Stata 系列推文:
作者:陈子厚 (中共广东省委党校)
邮箱:econometricalc@outlook.com
目录
计算莫兰指数的软件非常多,例如 Stata,详见连享会推文「Stata:面板数据的莫兰指数计算与散点图绘制-xtmoran」。本文主要借助 Matlab 实现莫兰指数的计算,有编程兴趣的同学,可以将此代码作为练习。
莫兰指数是全域莫兰指数 (Global Moran's I) 的简称,是澳大利亚统计学家帕特里克·阿尔弗雷德·皮尔斯·莫兰在 1950 年提出的,是最常用的空间自相关指标。其计算公式如下:
公式中
对于莫兰指数的检验,原假设
莫兰指数统计量的期望公式如下:
莫兰指数统计量的方差公式如下:
为了方便理解和实现 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
将数据另存为 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
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
Note:产生如下推文列表的 Stata 命令为:
lianxh 空间计量, m
安装最新版lianxh
命令:
ssc install lianxh, replace
免费公开课
最新课程-直播课
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
⛳ 课程主页
⛳ 课程主页
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会-常见问题解答:
✨ https://gitee.com/lianxh/Course/wikis
New!
lianxh
和songbl
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh