空间杜宾模型Matlab实现：Elhorst面板代码使用与报错排查指南

简介本资源是一套面向经济学、地理学及社会科学领域研究者的空间面板计量分析MATLAB工具包聚焦于空间杜宾模型SDM、空间滞后模型SAR与空间误差模型SEM的实现与调试特别适用于处理具有空间依赖性的面板数据并应对内生性问题。压缩包共57个文件主体为53个MATLAB函数.m涵盖权重矩阵特征值计算sar_eigs.m/sem_eigs.m、面板固定/随机效应估计sar_panel_FE.m/sem_panel_RE.m、LM检验与稳健标准误lmerror_robust_panel.m/lmlag_robust_panel.m、直接与间接效应分解direct_indirect_effects_estimates.m等核心模块另含3个WK1格式实证数据集如cigarette.wk1及1个MATLAB数据文件sp3.mat总大小仅143KB轻量实用。已有403人学习下载提供完整可运行代码框架、典型错误提示如panelcode error的上下文线索及多场景演示脚本demo*.m便于研究者快速复现模型、定位权重设定或估计收敛类问题并深入理解空间溢出效应的计量逻辑。 搞空间面板模型的人应该都见过Paul Elhorst那套Matlab代码。前阵子我从网上下了一个“New Elhorst Panel Code.zip”想跑空间杜宾模型Spatial Durbin Model也就是常说的SDM结果一上手就被panel code error折腾得够呛。今天不绕弯子把这套代码怎么用、跑杜宾模型时最容易踩的坑、以及各种报错怎么排查一次说清楚。如果你正准备用Matlab做空间面板估计或者已经在跟这套代码搏斗这篇文章应该能帮你省下好几个晚上。先说结论Elhorst这套代码本身写得非常规整绝大多数报错都不是代码的锅而是数据排序、权重矩阵、参数设置这三关没过。下面按我踩坑的顺序一层层拆开讲。1. 先搞清楚这套代码到底在算什么1.1 空间杜宾模型一个公式讲清楚空间杜宾模型的长相是这样的y ρWy Xβ WXθ ε其中y是被解释变量X是解释变量W是空间权重矩阵Wy叫空间滞后项WX叫解释变量的空间滞后。ρ是空间自回归系数β是X本身的系数θ是WX的系数ε是误差项。这个模型的厉害之处在于它同时考虑了三种空间关系本地区y受本地区X的影响β、受其他地区y的影响ρWy、还受其他地区X的影响WXθ。换句话说SDM既抓了空间溢出效应又抓了空间交互效应是空间计量实证里用得最广的模型之一。很多新手上来就问SAR空间滞后模型和SEM空间误差模型选哪个我的建议是如果没有强有力的理论依据直接用SDM做起点反而更稳。因为SAR只管WySEM只把空间相关塞进误差项而SDM把WX这一步补上了。漏掉WX相当于把解释变量的空间溢出效果硬生生忽略掉这在很多经济场景下是说不通的比如一个城市的环保政策会影响周边城市不只是影响自己的污染排放。不过注意SDM的估计不是普通OLS能搞定的因为存在Wy这一项OLS估计会不一致。Elhorst的代码用的是最大似然估计MLE这也是为什么他对数据排序和权重矩阵这么敏感——MLE里到处是矩阵运算一步错步步错。1.2 Elhorst面板代码包里有什么这套代码包里的核心文件大致分几类一是主估计函数比如sar_panel_FE.m、sem_panel_FE.m、sdm_panel_FE.m、sdm_panel_RE.m分别对应固定效应和随机效应下的不同模型二是检验函数比如lmsar_panel.m、lmsem_panel.m做LM检验lr_tests_panel.m做LR检验三是后估计函数比如direct_indirect_effects_2012.m专门做直接效应和间接效应分解。我得提醒一句网上流传的版本很多不同年份的包函数名和返回结构可能有细微差别。你下载的如果是“New Elhorst Panel Code”大概率是Elhorst本人2014年前后那版配套代码里面函数命名基本如上。但也有人解压后发现和另一份教程对不上号这时候别慌打开函数文件看注释以你手里这份的接口为准。这套代码适合谁适合已经有一份面板数据、想跑空间计量实证的硕士博士和研究人员。它不要求你会写MLE但要求你懂一点矩阵运算基础至少知道N×T和T×N的区别否则报错时你会完全懵。2. 跑代码前必须做好的三件事2.1 面板数据排序N×T还是T×N搞错全盘皆输这是整个代码包里最阴的一个坑没有之一。Elhorst的代码要求你的数据按“个体优先”的方式排列第一个个体的所有时期排完再排第二个个体以此类推。也就是先固定个体i再按时间t递增排列。举个例子你有3个地区、2个年份正确的排列顺序应该是地区1-2000年 地区1-2001年 地区2-2000年 地区2-2001年 地区3-2000年 地区3-2001年但很多人从Excel里整理数据时习惯按年份堆数据结果排成了“2000年三个地区、2001年三个地区”这在Stata里用xtset可能无所谓但在Elhorst这套Matlab代码里权重矩阵W是按个体顺序对应y的一旦数据顺序和W对不上估计出来的ρ直接就是错的而且错得没谱——正负号都可能反过来。我自己的血泪经验是拿到面板数据后先看id变量和year变量用sortrows按id、year排序然后再抽y和X。排序这件事花不了两分钟但能省下后面几小时的排查时间。我之前还遇过一个更隐蔽的情况数据里某些地区某年缺失别人用平均值填充导致同一地区相邻年份数据完全一样结果MLE迭代半天不收敛。遇到缺失值能删就删能用插值就用插值别用均值糊弄空间面板对数据质量的要求比普通面板更高。2.2 空间权重矩阵生成、标准化、以及必须避开的坑空间权重矩阵W是整个模型的心脏。它的维度必须是N×N这里的N是横截面单元数地区数不是NT。千万别把W写成NT×NT这是很多新手第一次跑代码必踩的雷。W的生成方式常见有三种邻接矩阵两个地区接壤记1否则记0、距离倒数矩阵距离越近权重越大、K近邻矩阵每个地区找最近的K个邻居。具体用哪一种取决于你的研究问题但在经济实证里最常见的是Queen邻接矩阵也就是只要两个地区有公共边界就算邻居。生成W之后必须要做的一步是行标准化也就是让每一行元素之和等于1。Elhorst代码包里自带normw函数调用方式很简单W normw(W_raw);行标准化的目的是让空间滞后项Wy可以被解释成“邻居的平均值”这样ρ的大小就能直接和空间依赖强度挂钩。没做行标准化之前ρ的含义会很模糊数值范围也不稳定。权重矩阵的坑相当多我单独列几个最常见的第一W矩阵对角线必须是0自己不能是自己的邻居第二如果某个地区没有任何邻居比如一个孤岛省份行标准化后会出NaN后面所有计算全崩第三W的个体顺序必须和y、X里个体的顺序一致很多人生成W时用的地区ID排序和主数据不一样结果跑出来结果异常还没察觉。我自己现在的工作流是用经纬度或者shp文件生成邻接关系后先画个图看一眼确认没有孤立点再检查一下sum(W,2)是否全为1最后才放心进模型。2.3 主程序参数设置固定效应、随机效应与info结构体Elhorst代码里估计函数通常长这样results sdm_panel_FE(y, X, W, T, info);其中y是被解释变量列向量长度NTX是解释变量矩阵NT×KW是N×N的行标准化权重矩阵T是时间期数info是一个结构体用来传各种参数。info里最重要的字段是info.flag它控制固定效应的类型。常见的约定是0表示无固定效应1表示空间固定效应个体固定效应2表示时间固定效应3表示空间和时间双向固定效应。具体到你自己下载的版本以代码注释为准但大体逻辑都是这个。为什么这个参数重要因为选错了固定效应设定模型结果的天差地别。空间固定效应相当于给每个地区一个自己的截距用来吸收那些不随时间变化但随地区变化的遗漏变量时间固定效应则吸收那些不随地区变化但随时间变化的遗漏变量。双向固定效应最严格但也最吃样本量。选择固定效应还是随机效应理论上应该做Hausman检验但很多实证文章直接按研究习惯来。我的建议是如果你的数据覆盖的个体几乎就是研究总体的全部比如全国31个省份直接用固定效应不用纠结随机效应如果你的个体是从总体中抽样得到、且希望推及总体才考虑随机效应。另外说一句info里还可能有个info.model之类的字段用来切换模型类型但即使有一般也不需要你手动改因为模型选择是靠调用不同函数实现的而不是靠改info。3. 杜宾模型估计与效应分解实操3.1 完整跑通一次SDM估计的步骤下面给一套可以直接照着改的完整流程假设数据在Excel里包含id、year、y、x1、x2、x3。clear; clc; % 第1步读取数据 data xlsread(panel_data.xlsx); % 假设第1列id第2列year第3列y第4到6列是x1、x2、x3 id data(:, 1); year data(:, 2); y data(:, 3); X data(:, 4:6); % 第2步按id、year排序确保个体优先排序 [~, idx] sortrows([id, year], [1, 2]); y y(idx); X X(idx, :); N length(unique(id)); % 个体数 T length(unique(year)); % 时期数 % 第3步生成并标准化空间权重矩阵 % 这里假设你已经有了原始邻接矩阵W_rawN×N W normw(W_raw); % 第4步设置info参数双向固定效应 info.flag 3; % 第5步估计空间杜宾模型固定效应 results sdm_panel_FE(y, X, W, T, info); % 第6步查看结果 results.beta results.teta results.rho跑完之后results里会给你几个关键输出results.beta是X的系数results.teta是WX的系数results.rho是空间自回归系数。这三个数怎么看β看直接解释变量本身的影响θ看解释变量从邻居传导过来的影响ρ看被解释变量在地区间的依赖程度。ρ显著为正说明存在正向空间溢出ρ为负说明地区间是竞争或挤出关系。有一点要提醒SDM结果不能只报回归系数正文里一定要给直接效应和间接效应分解因为SDM的系数本身不直接等于边际效应。这是审稿人最爱挑的点下面细说。3.2 直接效应、间接效应与总效应别只盯着回归系数在空间杜宾模型里某个解释变量x的变化影响的不只是本地区的y还会通过空间传导影响其他地区的y。要量化这种影响得用LeSage和Pace提出的偏微分方法把总效应拆成直接效应和间接效应。理解方式很简单直接效应是“x变了对本地区y的平均影响”这里面已经包含了空间反馈比如本地区x提高带动邻居y提高邻居y又反馈回来影响本地区y间接效应是“所有其他地区x同时变对本地区y的平均影响”也就是真正的空间溢出。Elhorst代码包里提供了现成的分解函数用法大致是direct_indirect_effects_2012(results, W, N, T);它会输出每个解释变量的直接效应、间接效应和总效应以及对应的t统计量。记得在论文里报告标准差或p值光给点估计会被质疑。我见过不少同学跑完SDM把results.beta里的系数当成边际效应直接写进论文这是不对的。举个具体例子某研究里x的β是0.3看着不大但直接效应算出来是0.42间接效应是0.18总效应0.6。如果你只看β会严重低估政策变量的总影响。效应的数量级和显著性跟β是有差异的必须单独报告。3.3 模型选择LM检验、LR检验与Hausman检验怎么配合不要拿到数据就直接说“我跑SDM”应该有一串检验流程。常用的三件套是LM检验、LR检验、Hausman检验。第一步用LM检验判断空间相关主要来自哪个方向。Elhorst代码里的lmsar_panel.m检验“空间滞后”方向lmsem_panel.m检验“空间误差”方向。这两个检验都有稳健版本robust LM报告时最好同时给普通版和稳健版。第二步用LR检验判断SDM能不能简化成SAR或SEM。逻辑是这样的如果SDM里所有θWX的系数联合为0那SDM就退化成SAR如果θ ρβ 0这个约束成立那SDM退化成SEM。Elhorst包里lr_tests_panel.m能同时给出这两个检验的结果分别对应“SDM能否简化成SAR”和“SDM能否简化成SEM”。如果两个原假设都被拒绝就用SDM如果只有一个被拒绝可以选更简约的模型。第三步才是Hausman检验判断固定效应还是随机效应。这块很多软件都能跑Matlab里也可以用现成命令或手写。我个人的实操顺序是先跑一个不带空间项的普通面板看LM检验结果初步判断空间依赖形态然后直接跑SDM做LR检验确认是不是必须用SDM最后用Hausman定效应类型。这套流程走完模型选择环节基本无懈可击。4. panel code error高频报错与排查实录4.1 维度不匹配出现频率最高的报错如果你只记一条记住这条报错里只要出现“Matrix dimensions must agree”九成是维度问题而且八成出在W上面。具体来说检查三件事一是W是不是N×N拿size(W)看一眼如果发现W是NT×NT或别的尺寸赶紧改二是y是不是NT×1、X是不是NT×K很多人用xlsread读数据时会带出表头行导致y和X长度差一格三是T和N的赋值对不对T是时间长度不是年份最大值N是个体数不是数据总行数。还有一次我帮一个师弟排查代码里写的是N length(data)结果W是31×31y是310×1N被设成了310后面所有关于N的运算全乱套。这种错误很好修但排查起来很费时间因为报错不一定直接指向N而是指向后面某个矩阵乘法。4.2 特征值与矩阵运算错误MLE里经常要对某个包含W的矩阵做特征值分解所以报错里出现eig、eigenvalue、complex相关字样时别慌优先查W。最常见的一个问题是W没有标准化或者W里有NaN。一旦某一行全为0孤立点normw之后会出现0/0得到NaN后面的特征值计算直接崩。解决办法是回到权重矩阵生成那一步检查每个地区是否至少有一个邻居。实在有孤立点考虑改用距离权重把远处地区也连上或者直接删掉那个样本。另一个不太容易注意的问题是W不对称。邻接矩阵本身应该是对称的但如果你手工构造距离权重很容易造出一个不对称的W。虽然行标准化后每行和为1不一定报错但特征值会出现复数后续计算会出现奇怪的虚部。排查办法在跑模型前加一行assert(isequal(W_raw, W_raw))提前把不对称问题暴露出来。4.3 权重矩阵引发的诡异结果有些时候代码不报错但结果看着就不对这种情况最让人头大。我归纳过几种典型症状第一种ρ无比接近1或者大于1。这通常意味着空间依赖被过度放大常见原因就是W和数据顺序不对应或者W的行标准化没做对导致Wy这一项几乎等于y本身。第二种所有系数都巨大符号还反直觉。这种情况往往是个别地区是异常点而且异常点恰好又有很高的权重连接。可以画一下残差图看看有没有哪个地区残差特别大。曾经有个案例某个城市的增长率极端高是其他城市的几十倍加进去之后所有系数全乱删掉后模型就正常了。第三种结果对W的生成方式特别敏感。换一种权重矩阵设定显著性跟着变。这不是代码bug而是模型本身对W依赖太强稳健性欠佳。遇到这种情况我习惯于做敏感性分析邻接矩阵、距离倒数矩阵、K近邻矩阵各跑一遍如果核心结论不变再动笔写论文。4.4 结果异常不收敛、系数爆炸、符号不对最后说说MLE不收敛的事。Elhorst代码用的迭代优化有时候会报“Iteration limit exceeded”或者结果里出现NaN。先别急着改代码按顺序排查数据有没有标准化当X的量级差异特别大比如一个变量是0到1的比率另一个是几百万的人口数矩阵条件数会非常差迭代很难收敛。解决办法是把X做归一化或者取对数尤其注意别把几个数量级差太多的变量直接放一起。再检查一下T和N的设定是否跟真实数据一致。曾经有人数据里有个地区只有8年其他地区10年他硬把T设成10结果缺失时段的NaN直接传染给后面所有运算。处理不平衡面板时要么补齐数据要么换成专门处理非平衡面板的版本。还有一类情况是符号整体反了。如果所有系数符号都跟预期相反很可能是权重矩阵的方向反了或者数据里y/X的顺序没对齐。这几个问题不报错但结果属于“错得很有规律”需要靠你对业务的判断力来发现。我的一般做法是先用一个小模拟数据或者只用前两年数据快速验证一遍程序能跑通、符号能对上再上全样本。这套自检习惯帮我避免了好几次大返工。再说一个冷门但真实存在的坑有的zip下下来解压后是空文件夹或者里面少了几个关键函数文件比如normw.m没被包含进去。如果运行报“Undefined function normw”或者“Undefined function sdm_panel_FE”优先检查是不是文件缺失而不是代码写错。去Elhorst的页面把完整包重新下一遍或者从代码包里另外一个文件里把normw函数单独抠出来都能解决。根据我个人的使用经验这套代码在数据规范、权重矩阵正确的前提下跑起来是很稳定的。真正花时间的从来不是估计那一下而是前期的数据准备和后期对效应的解读。如果你正准备做空间杜宾模型我建议把上面这三件事先做扎实数据排序、权重矩阵标准化、效应分解报告。这三关过了后面基本一马平川。本文还有配套的精品资源点击获取