简介:本资源是一套面向经济学、地理学及社会科学领域研究者的MATLAB空间计量建模工具包,聚焦面板数据下的空间杜宾模型(SDM)实现与调试,特别适用于需处理空间依赖性与内生性问题的实证分析场景。压缩包共57个文件,主体为53个MATLAB函数(.m),涵盖空间滞后模型(SAR)、空间误差模型(SEM)及空间杜宾模型(SDM)的面板估计、似然计算、直接/间接效应分解、LM检验、稳健标准误计算等核心功能;另含3个WK1格式的空间权重矩阵文件与1个MAT文件,支撑实证数据加载与权重构建。目前已有403人学习下载,资源结构清晰、模块化程度高,提供从数据预处理(如demean.m)、模型拟合(如sar_panel_FE.m、sem_panel_RE.m)、诊断检验(如panel_test.m、LMsarsem_panel.m)到结果解读(如direct_indirect_effects_estimates.m)的完整技术链,附带多个演示脚本(demo*.m)与实证案例(cigarette.wk1等),是掌握Elhorst空间面板方法的实用代码参考。 几个月前我收到一条读者私信,内容特别简短:下载了New Elhorst Panel Code.zip,准备跑空间杜宾模型(SDM),结果panelcode一运行就报error,卡了三天没进展。这个场景我太熟了。Elhorst的空间面板数据代码在空间计量圈子里流传极广,几乎所有做空间杜宾模型、空间滞后模型、空间误差模型的实证研究者,手里都有一份从某个学术群或网盘里转存来的"New Elhorst Panel Code.zip"。但真正能第一次就跑通的人,十个里面未必有一个。多数人不是倒在模型原理上,而是倒在最基础的代码调用、数据组织和路径配置上。
这篇文章就围绕Elhorst这套面板代码,把空间杜宾模型的MATLAB实现从头到尾捋一遍。重点解决两件事:一是怎么正确地把代码跑起来,二是跑的时候那个"panelcode error"到底是怎么来的、怎么查、怎么修。全程按我实际操作的思路来写,涉及的具体函数名、脚本写法、排查命令,都是我长期用下来比较稳定的方案,你可以直接照抄。
1. 先把这套代码和杜宾模型的关系搞清楚
1.1 Elhorst代码包到底装了什么
Paul Elhorst是荷兰格罗宁根大学的空间计量经济学教授,他在个人主页上长期维护一套空间面板数据模型的MATLAB程序。网上流传的"New Elhorst Panel Code.zip",主要就是他主页提供的那套更新版面板代码。为什么这套代码在学术界传播度这么高?因为它把LeSage和Pace在《Introduction to Spatial Econometrics》里推导的极大似然估计、空间效应分解等理论,变成了能直接运行的MATLAB函数。你不需要自己写似然函数的迭代代码,只要准备好面板数据y、解释变量X、空间权重矩阵W,就能估计出SAR、SEM、SDM、SAC等一系列空间面板模型。
代码包内部结构因版本而异,但核心文件一般包括:
- f_sarpanel.m、f_sempanel.m、f_sdmpanel.m、f_sacpanel.m这些模型估计主函数
- f_xxx_direct_indirect_effects.m这类效应分解函数
- 一个或几个demo脚本,告诉你如何组织数据并调用函数
有点讽刺的是,新版代码包里demo脚本往往是"残缺"的。它有函数原型、有参数说明,但缺少一个完整的main脚本示例,需要你自己把数据整理成特定格式再调用。这一步恰恰是新手最容易卡住的地方:代码本身没错,但你的数据格式不符合函数签名要求,于是MATLAB甩出一串红色报错。
1.2 杜宾模型究竟在做什么
空间杜宾模型(Spatial Durbin Model,SDM)在普通线性回归基础上加了两样东西:被解释变量的空间滞后项,以及解释变量的空间滞后项。用公式表示就是:
y = rho * W * y + X * beta + W * X * theta + epsilon
其中rho是空间自回归参数,W是空间权重矩阵,beta是解释变量系数,theta是解释变量空间滞后项的系数。
这个模型的优势在于,它把"邻居的影响"拆成两条路径。一条路径是邻居的y变动直接带动本地y变动,比如某个城市房价上涨,周边城市房价跟着被拉高;另一条路径是邻居的x变动通过溢出效应影响本地y,比如周边城市人均收入提高后,通过跨区域消费和投资带动本地房价。如果你只关心"有没有溢出效应"以及"溢出效应的规模有多大",SDM是比SAR更稳妥的起点:SAR只允许y的空间溢出,SEM则把空间相关性塞进误差项,都不如SDM解释得直接。
但这引出一个关键操作要求:SDM估计出的beta和theta不能直接当作边际效应来报告。必须把总效应拆成直接效应(本区域x对本区域y的影响)和间接效应(本区域x对其他区域y的平均影响)。Elhorst代码里专门有函数做这个分解,等到第3节我具体讲怎么用。
1.3 panelcode error 这个关键词从哪来
搜索词里出现"panelcode error",其实"panelcode"并不是某个官方函数名,而是大家在论坛、网盘、学术交流群里对"面板数据代码"的简称。当一个人下载了"New Elhorst Panel Code.zip",运行demo时MATLAB命令窗口报出一条错误,他不知道该搜什么,就用"panelcode error"去找答案。这个检索词背后,反映的是比单条报错更深一层的问题:多数人拿到这套代码后,不知道代码的运行机制是什么,也不知道该怎么查错。埃尔霍斯特的代码并不是成熟的商业软件,它是一套研究用程序,需要你理解数据组织方式、函数调用约定和估计原理。任何一处偏差都可能让代码停在半路。
2. 运行前的准备工作:环境、数据、权重矩阵
2.1 MATLAB环境与依赖工具箱检查
在解压zip包之前,先确认MATLAB版本。Elhorst新版面板代码基于统计工具箱和优化工具箱,建议R2016a以上,我在R2020a、R2022b上都实测过,可以正常运行。如果你还停留在R2014a这种老版本,部分函数语法可能不兼容,比如某些矩阵运算的写法在新旧版本中存在差异。
更关键的是jplv7工具箱。Elhorst很多函数内部调用了LeSage的空间计量工具箱jplv7中的子函数,比如norm_rnd、stdn等。如果你只下载了Elhorst的zip包,没有把jplv7加入MATLAB搜索路径,那么运行时马上会报"未定义函数或变量"。
验证环境是否就绪,最简单的方式是在命令窗口输入:
which f_sdmpanel which norm_rnd如果两行命令都返回路径,说明环境和依赖都正常。如果返回"未找到"或空白,就要添加路径:
addpath(genpath('你的jplv7文件夹路径')); addpath(genpath('你的Elhorst代码文件夹路径')); savepath;savepath一定要执行,否则下次启动MATLAB又要重新添加。
2.2 数据格式的硬性要求
Elhorst代码对面板数据顺序有极其严格的要求。假设有N个截面单元(比如30个省份),T个时期(比如2010到2022年,共13年),那么:
- y必须是N*T行、1列的列向量,排列顺序是:第1个截面单元的T个时期数据排在最前,接着是第2个截面单元的T个时期数据,依此类推
- x是N*T行、K列的矩阵,每一列对应一个解释变量,排布顺序与y一致
- W是N行N列的空间权重矩阵,与时间维度无关
这个"先截面后时间"的排序方式,与大家日常习惯的"先时间后截面"(也就是每年报表里各省份按行排、年份按顺序堆叠)完全不同。很多人把数据按Excel里的原始顺序导入,不重排,代码运行一般不报错,但估计结果完全错乱。这是最隐蔽的坑:不是报文化错,而是结果错。
快速检查排序是否正确,可以用这样一段代码:
T = 13; % 时期数 N = 30; % 截面数 disp(y(1:T)); % 应该是第一个截面单元的T个时期观测 disp(y(T+1:2*T)); % 应该是第二个截面单元的T个时期观测如果y的前T个数据并非第一个个体按期排列,你就需要重新整理数据顺序。最稳妥的办法是生成一列"个体编号"和一列"年份",在Excel里按个体编号升序、年份升序排序后再导出。
2.3 空间权重矩阵W的构造与标准化
W是空间计量的灵魂,也是出错重灾区。Elhorst代码要求空间权重矩阵必须满足三个条件:
- N行N列方阵
- 对角线全为0
- 每行之和等于1,即行标准化
常见W构造方法有三种。第一种是邻接矩阵:两个区域有共同边界则记为1,否则为0,然后行标准化。第二种是距离矩阵:用经纬度计算球面距离或欧氏距离,取距离倒数或距离倒数的幂。第三种是经济距离矩阵:基于人均GDP等经济指标差异来定义"经济上的邻居"。
无论哪种方法,最后都要做行标准化。下面是一个典型的构造流程:
% 假设W_raw是原始构造的N*N邻接矩阵,对角线为0 W = W_raw ./ sum(W_raw, 2); % 检查是否标准化成功 disp(sum(W, 2));如果输出结果中有NaN或0,说明W_raw中有一行全为0。这种情况常见于邻接矩阵:某个地区没有任何相邻地区,比如海岛城市或偏远省份。此时行标准化时分母为0,产生NaN。解决方法是改用"最近K个邻居"矩阵,或者直接使用距离权重矩阵。
我在处理地级市数据时频繁遇到这个问题。一个城市如果在地理上孤立,就会让整个N*N矩阵出现坏行,进而导致极大似然估计无法计算。
3. 杜宾模型代码实操:从main脚本到结果解读
3.1 最小的main脚本长什么样
把数据处理成Elhorst要求的格式后,main脚本其实很短。以SDM面板模型为例:
% 清理环境 clear; clc; % 加载数据 % 假设y为N*T x 1列向量,x为N*T x K矩阵,W为N x N行标准化权重矩阵 load('mydata.mat'); N = 30; % 截面数 T = 13; % 时期数 % 调用SDM面板估计 result = f_sdmpanel(y, x, W, N, T); % 输出直接效应与间接效应 f_sdm_direct_indirect_effects(result); % 显示主要结果 disp(result.beta); disp(result.rho); disp(result.tstat);这里有一个需要留意的点:f_sdmpanel这个函数名在不同版本的Elhorst代码里可能不同。旧版可能叫f_sdm,新版可能叫f_sdmpanel,还有一些版本带了后缀。打开你下载的代码包,看看到底存在哪个.m文件,以实际文件名为准。如果不确定,在命令窗口输入dir *.m,列出所有函数文件。
3.2 参数设置的细节
f_sdmpanel的函数签名一般是:
function results = f_sdmpanel(y, x, W, N, T, info)info是一个可选的优化参数结构体,常用字段有三个:
- info.lflag:似然函数逼近方式。lflag=0用精确似然,小样本、中等样本推荐;lflag=1用近似似然,大样本时速度快很多
- info.rmin、info.rmax:空间自回归参数rho的搜索区间。默认通常是(-1, 1),有时可以放宽到(-1.5, 1.5)
- info.convg:收敛精度,默认1e-5,一般不需要改动
对于常规省级面板,我建议直接用:
info.lflag = 0; info.rmin = -1; info.rmax = 1; result = f_sdmpanel(y, x, W, N, T, info);如果数据量特别大,比如N=300、T=20,那么NT为6000行,精确似然计算涉及反复计算NN矩阵的行列式,速度会很慢。这时把lflag设为1,能明显提速。不过要注意,近似似然的结果与精确似然可能略有差异,论文里建议主结果用精确似然,稳健性检验里可以报告近似似然的结果。
3.3 模型结果的字段与效应分解
Elhorst代码运行结束后,result是一个结构体,常用字段包括:
- result.beta:解释变量系数估计
- result.rho:空间自回归参数
- result.tstat:各系数的t统计量
- result.yhat、result.resid:拟合值和残差
- result.lik:极大似然函数值
- result.rsqr:拟合优度
但实证报告中不能只报beta。SDM的核心结论来自直接效应和间接效应。直接效应反映本区域解释变量变化对本区域被解释变量的平均影响,间接效应反映本区域解释变量变化对其他区域被解释变量的平均影响,也就是空间溢出效应。Elhorst在代码包里提供了专门的效应分解函数:
% 直接调用即可 f_sdm_direct_indirect_effects(result);这个函数会在命令窗口输出各变量的直接效应、间接效应和总效应,以及对应的t统计量。写论文时,把这些效应值整理成表格,比单纯列beta更有说服力。审稿人看到你没有做效应分解,一眼就能看出你对SDM的理解还停留在表面。
4. panelcode error:高频报错与排查实录
4.1 路径缺失与依赖函数找不到
报错信息通常是"未定义函数或变量 'norm_rnd'",或者"未定义函数或变量 'f_sdmpanel'"。
原因很简单:jplv7或Elhorst代码目录没有加入MATLAB路径。这种问题占初学者报错的三分之一以上。
解决方法:
addpath(genpath('D:\work\jplv7')); addpath(genpath('D:\work\elhorst_panel')); savepath;注意genpath会自动添加文件夹下所有子文件夹,避免漏掉jplv7里嵌套的子目录。
4.2 矩阵维度不一致
报错信息可能是"矩阵维度必须一致",也可能是"索引超出数组边界"。
原因几乎都是数据格式问题:y的行数不等于NT,或者W的维度不是NN。
一个特别常见的低级错误是把W扩展成了NT行、NT列的矩阵。W表示的是截面单元之间不随时间改变的空间关系,它永远是N*N。无论你的面板有多少期,W都只有一张。
排查方法:
disp(size(y)); % 应该输出 N*T, 1 disp(size(x)); % 应该输出 N*T, K disp(size(W)); % 应该输出 N, N disp(N * T); % 必须等于 size(y,1)只要这四行输出对不上,就一口气顺着数回去,重排数据。
4.3 NaN和Inf混进了数据
报错信息:"矩阵包含 NaN 或 Inf"。
原因有两类:一是数据源本身存在缺失值,没有清洗就直接进模型;二是空间权重矩阵行标准化时出现某一行全为0,导致NaN。
处理方式:
- 缺失值尽量用插值补齐,或者干脆剔除该截面单元,不要留NaN
- 权重矩阵某一行为0,需要修改W构造方法。比如改用距离权重矩阵,或者用"最近K个邻居"法生成邻接关系
我在处理中国地级市数据时,发现"孤岛"问题特别常见。一个海岛城市按地理邻接可能找不到任何邻居,行标准化直接产生NaN。换成距离矩阵后,每个城市都能算出与最近城市之间的距离,问题自然解决。
4.4 极大似然迭代不收敛或rho跑边界
这种问题通常不报红色error,而是出现警告,或者估计结果里rho非常接近1(比如0.9999),beta的符号明显异常。
排查方向:
- 检查W有没有完成行标准化
- 检查X和WX之间是否存在严重共线性。SDM同时包含X和WX,如果某个变量与其空间滞后项的相关性非常高,比如超过0.95,估计就会不稳定
- 尝试把模型退化为SAR,去掉WX项,看结果是否稳定
如果rho跑到边界,一个可行的做法是修改info.rmin和info.rmax范围,比如放宽到(-1.5, 1.5),给优化器更大的搜索空间。但如果放宽后rho依然贴着边界,就要怀疑数据或模型设定本身有问题。
4.5 计算慢得像死机
空间面板最大似然估计的复杂度很高,每一步迭代都要处理N*N矩阵的行列式。N=300时,单次迭代可能就要数十秒,T又较大时,整个过程跑几个小时很常见。很多初学者以为代码卡死了,其实它只是慢。
优化建议:
- 把info.lflag设为1,使用近似似然
- 用稀疏矩阵存储W
- 尽量把数据放在内存里,不要从Excel实时读取
- 用MATLAB的batch模式在后台跑,别一直盯着进度
4.6 一个快速定位错误的通用方法
如果报错信息看不懂,或者不知道错在哪一行,有个非常实用的技巧:在main脚本开头加一行调试命令。
dbstop if error加上这行后,程序一旦报错,MATLAB会自动停在出错的那一行,进入debug模式。在工作区里直接查看所有变量的尺寸、值、是否为NaN,比反复读错误信息直观得多。我用这个方法帮很多学生排查Elhorst代码,90%的问题都在一分钟内定位:不是y的行数不对,就是W没有标准化,或者某个数据列混入了文本导致矩阵变了类型。
5. 避坑心得与扩展建议
5.1 我实际踩过的几个特殊坑
第一,数据文件从Excel读入时,readmatrix可能把表头当成了数据。一个典型的特征是MATLAB提示"数据必须是数值型",而你的Excel第一行恰好是中文列名。建议导出数据时删掉表头行,或者用readmatrix指定Range。
第二,如果你使用中国省级面板数据,注意省级行政区划代码在不同年份可能存在调整。比如某些年份县级市升级为地级市、省直辖县变动等,这类变化会让两个年度的截面单元口径不一致,导致合并后的面板数据排序错乱。最稳妥的做法是:先构建一个稳定的省份名单,再按名单顺序整理每一年的数据,最后堆叠。
第三,Elhorst代码包在多次转存后,可能出现同名函数的新旧版本并存。比如文件夹里同时有f_sdmpanel.m和f_sdmpanel_new.m。运行demo前先确认demo调用的是哪个版本,避免新旧版本混用,导致结果异常。
5.2 跑通之后还要做什么
跑通SDM只是实证的第一步。一篇能发表的论文还需要做三件事:
- 稳健性检验:换空间权重矩阵(邻接换成距离)、换模型(SDM退化成SAR和SEM)、换解释变量定义
- 效应可视化:把直接效应和间接效应画成柱状图或地图
- 与其他模型比较:用LR检验或AIC/BIC判断SDM是否显著优于SAR和SEM
Elhorst代码主要解决模型估计,可视化和模型比较部分需要你自己写。我习惯把效应分解的结果导出到Excel,再用Python的matplotlib或ArcGIS做空间分布图,这样既保留MATLAB的估计精度,又能让图更灵活。
5.3 如果不一定用MATLAB
现在Python的空间计量生态逐渐成熟,libpysp、spreg、spmodel这些库已经能处理多种空间面板模型。如果你的主要需求不是复现特定文献,而是开启新的实证项目,从Python入手可能更省力。但Elhorst代码在学术界的影响力和被引用频率仍然很高,很多时候你下载的参考代码、课上演示、师兄师姐的模板都是基于这套MATLAB代码。在这种情况下,掌握它仍然有不可替代的价值。
最后分享一点个人经验
我当年第一次用Elhorst代码跑SDM,整整折腾了一周。最后发现问题极其简单:数据排序是按时间而不是按截面排的,W完全对不上。那段经历让我养成了两个习惯。第一个习惯是:拿到任何空间计量代码,第一件事不是跑模型,而是先用size、sum、max这些命令把y、x、W的结构彻底检查一遍。第二个习惯是:遇到报错不急着改代码,先看数据结构,再看依赖路径,最后才追溯到模型设定。空间计量代码跟普通回归代码最大的区别,就是多了一个空间权重矩阵W,而W不直接出现在回归方程里,它像一个隐形的第三方,一旦出错,代码的报错信息往往指向莫名其妙的行列式或矩阵运算,让你摸不着头脑。希望这篇文章能让你少走几天弯路。
本文还有配套的精品资源,点击获取