★ FEATURED ARTICLE 简介面向经济地理与区域科学领域的研究者这份Elhorst空间面板数据模型新版本代码包可用来处理区域间的空间自相关性与时间序列动态性大幅简化模型设定与参数估计流程。压缩包共含三十二个文件以二十七个m脚本为主覆盖空间面板模型的主要估计函数、直接与间接效应分解、空间滞后/误差项设定以及LM检验等核心模块另附三个wk1示例数据文件和两个asv备份文件整个包体仅九十八KB轻量易用。已有六十七人学习使用适合需要开展空间计量分析的高年级学生、科研人员与政策分析师。代码包整合了Elhorst教授的新面板代码配合Cigarette等经典数据集可在MATLAB中直接运行复现完整的空间面板模型分析流程帮助使用者理解空间权重矩阵构建、结果解读与诊断检验等关键技术。适用于政策空间效应评估、产业聚集机制识别、污染物扩散模式分析等应用场景是经济地理学者与区域分析师的实用工具。1. 翻遍论坛才凑齐的 elhorst_model_new空间面板模型就从这套 Matlab 代码跑起做空间计量的人迟早会撞上 Elhorst 这个名字。elhorst_model_new 这个 rar 包是我当年从某学术论坛的旧帖子里翻出来的解压下来十几个 .m 文件没有说明书只有一个 readme 用半页纸列了函数名。这套代码覆盖了空间面板模型最常用的估计路径SAR、SEM、SDM、SAC以及对应的固定效应、随机效应版本估计核心是极大似然ML。如果你手里有一批 N×T 的面板数据需要同时处理空间依赖和个体异质性这个包就是最省事的起点。新手拿它入门跑通第一个空间杜宾模型很合适老手也常留一份用来交叉验证论文结果。这篇笔记按我自己拆包的使用顺序写先讲结构再讲跑通最后讲参数与踩坑。2. 拆开 rar 看内部结构这些 .m 函数各自负责哪一类模型2.1 文件分类与命名规律解压 elhorst_model_new 之后你会看到一堆命名高度规律化的 .m 文件。只要看懂命名规则整个包的结构就一目了然。函数名一般是三段式模型缩写 数据形式 效应类型。我把最常用的一批函数整理成了下面这张表方便你对照手头的压缩包逐项确认。函数名模型类型效应形式返回值要点sar_panel_FE.m空间自回归 SAR固定效应beta、tstat、rho、liksar_panel_RE.m空间自回归 SAR随机效应在 FE 基础上增加 etarhosem_panel_FE.m空间误差 SEM固定效应beta、tstat、lambda、liksem_panel_RE.m空间误差 SEM随机效应同上sdm_panel_FE.m空间杜宾 SDM固定效应含直接/间接效应字段sdm_panel_RE.m空间杜宾 SDM随机效应含直接/间接效应字段sac_panel_FE.m空间自回归 误差组合 SAC固定效应包含 rho 和 lambdalmlag_panel.m空间滞后 LM 检验无LM 统计量、p 值lmerror_panel.m空间误差 LM 检验无LM 统计量、p 值lratios.mLR 似然比检验无LR 值、自由度waldtests.mWald 检验无Wald 值、p 值normw.m权重矩阵行标准化无标准化后的 W命名规律很直白sar 对应空间滞后项sem 对应空间误差项sdm 同时包含滞后项和误差项的解释变量交叉项sac 是滞后项和误差项同时存在的混合模型。面板数据专用后缀是 panel固定效应是 FE随机效应是 RE。如果你解压出来的文件和我列的不完全一样也不用慌。不同渠道流传的版本会在检验函数上做增删但核心估计函数基本都在。先把表里的 sar_panel_FE、sdm_panel_FE、normw 找齐这套代码就能用起来。2.2 ML 估计的核心逻辑为什么不能用普通 OLS 硬跑空间面板模型不能用 OLS 直接估计原因是空间滞后项 Wy 和误差项里的 Wε 都包含了内生性的空间反馈。以 SAR 模型为例y ρWy Xβ ε把 Wy 移到等式左边再求逆任意一个地区的扰动都会通过空间传导影响所有其他地区再反作用于该地区自身。这种反馈 loop 导致 OLS 估计量失真必须用极大似然方法处理。Elhorst 这套代码的 ML 估计过程分三步。第一步对原始数据进行 within 变换demean把个体效应 α 从方程里消掉避免直接估计 N 个虚拟变量导致的维度爆炸。第二步对 ρ 在可行区间内做网格搜索每给一个 ρ就用广义最小二乘求出 β 和 σ² 的条件估计对应一个对数似然值。第三步在所有候选 ρ 里挑对数似然最大的那个组合再用牛顿迭代精修。代码里这套步骤是封装在黑匣子里的你不需要自己写似然函数。但理解这个流程很重要因为后面所有参数的解读都建立在 ML 估计的基础上rho 代表空间依赖强度sige 是误差方差lik 是收敛后的对数似然值模型比较时经常要用到它。2.3 为什么选这一版而不是其他空间计量工具包空间计量 Matlab 工具包不止 Elhorst 这一套常见的还有 LeSage 的 jplv7 和后来维护的 spreg。我的习惯是论文复现和教学场景优先用 Elhorst原因有三。第一jplv7 的面板函数比如 sar_panel偏向早期实现固定效应靠手动构造虚拟变量数据量稍大就跑得吃力Elhorst 这版把 demean 过程内置N 到几千都不会内存爆炸。第二Elhorst 2014 版之后的代码在 SDM 模型里直接返回直接效应、间接效应和总效应的点估计与 t 值而 jplv7 需要自己写 delta method 去推标准误。第三这套代码的 LM、LR、Wald 检验函数齐全做模型选择时不用在多个工具包之间来回切换。当然这不代表 Elhorst 全方位优于其他包。LeSage 的代码在截面空间模型的贝叶斯估计上更强如果你要做 MCMC 抽样还是得回到 jplv7。工具选型本来就没有银弹我的原则是面板数据 ML 估计 需要效应分解时先拿 Elhorst 跑一版基线结果。3. 跑通第一个模型数据格式、W 矩阵构建与函数调用的完整流程3.1 数据怎么摆y、x、W 三者缺一不可Elhorst 代码对输入格式有严格要求这一步错了后面全是白费。先说结论y 是 NT×1 的列向量x 是 NT×k 的矩阵k 是解释变量个数不含常数项W 是 N×N 的空间权重矩阵。排序规则是代码里最容易翻车的地方。包内的估计函数默认面板数据按“个体为块、个体内部按时间升序”排列——前 T 行属于个体 1接着 T 行属于个体 2以此类推。如果你的原始数据是按时间整理的前 N 行是第一期所有个体必须先做一次排序再喂给函数。% 原始数据: id(个体编号), year(年份), y(因变量), x1, x2(解释变量) % 目标排序: 个体1全部年份 - 个体2全部年份 - ... % 先按个体分组, 组内按年份升序 data sortrows(data, [1 2]); % 第1列是id, 第2列是year % 提取变量 y data.y; % NT×1 列向量 x [data.x1, data.x2]; % NT×2 矩阵, 解释变量 N length(unique(data.id)); % 个体数 T length(unique(data.year)); % 期数 % 检查维度 assert(length(y) N*T, y长度必须等于N*T);排序是最不起眼却最容易出问题的操作。很多情况下估计程序不报错、迭代也正常收敛但因为 y 和 x 的个体对应关系错位估计出的系数符号和大小完全偏离理论预期。我一般在排序后专门抽几个个体核对一次确认每个个体的 T 行数据完整且连贯。数据缺省值的处理也要提前做好。面板数据最常见的坑是某个个体某期缺失直接删行会破坏排序结构。先做插补或者把对应的个体整个剔除务必保证每个个体的期数严格相同否则 reshape 时会报维度错误。3.2 跑通 SAR 固定效应模型数据准备好之后第一个模型我建议从 SAR 固定效应开始。SAR 只包含一个空间滞后项结构最简单跑通它就能验证数据格式和 W 矩阵都没问题。% 载入权重矩阵 W, W 是 N×N 稀疏矩阵, 对角线为0 % 常见做法是从邻接表构造, 或者读入外部生成的矩阵文件 load(W.mat); % 行标准化: 每行除以行和, 保证每行加总为1 W normw(W); % 调用固定效应 SAR 估计 results sar_panel_FE(y, x, W, T); % 查看核心输出 fprintf(rho %.4f\n, results.rho); fprintf(beta1 %.4f (t %.4f)\n, results.beta(1), results.tstat(1)); fprintf(beta2 %.4f (t %.4f)\n, results.beta(2), results.tstat(2)); fprintf(log-likelihood %.4f\n, results.lik);normw 函数做的是行标准化标准化后的 W 每一行和为 1这样空间滞后项 Wy 的系数 ρ 可以解释为邻接地区的加权平均影响。如果你的原始 W 是对称的邻接矩阵行标准化会破坏对称性这是刻意为之属于空间计量里的标准处理方式。调用函数时传入的 T 必须和数据实际的期数一致。函数内部会拿 T 把 y 重新排列成 T×N 的矩阵然后做 demean。如果 T 传错了内部维度检查和排序会全部错乱最典型的症状是报错“Dimensions must agree”或者估计结果出现大面积的 NaN。跑通后重点关注三个数rho 是否在 0 到 1 之间、beta 的符号是否和理论预期一致、log-likelihood 是否是一个合理的负数。这三个数任一出现异常先回数据格式和 W 矩阵排查。3.3 从 SAR 到 SDM自动生成 WX 与效应分解SDM空间杜宾模型在 SAR 的基础上加入了 WX 项即邻居的解释变量也会影响本地区因变量。模型表达式为 y ρWy Xβ WXθ ε。Elhorst 代码的方便之处在于你不需要手动构造 WX函数内部会基于 W 和 x 自动生成并返回对应的系数向量。% 调用固定效应 SDM 估计 results_sdm sdm_panel_FE(y, x, W, T); % 输出系数: 第一个是beta, 后面是WX的theta fprintf(rho %.4f\n, results_sdm.rho); fprintf(beta ); disp(results_sdm.beta); fprintf(theta (WX系数) ); disp(results_sdm.theta); % 关键: 直接效应、间接效应、总效应 fprintf(直接效应:\n); disp(results_sdm.direct); fprintf(间接效应:\n); disp(results_sdm.indirect); fprintf(总效应:\n); disp(results_sdm.total);这里有个极易混淆的点我见过不少用户手动把 WX 加进 x 矩阵再跑 SAR结果 SDM 的 theta 始终不收敛。这是典型的重复构造错误——SDM 函数内部会自己生成 WX如果外部再喂一份等于模型里同时出现了两组 WX设计矩阵发生共线性估计自然失败。正确做法是x 里只放原始解释变量WX 交给函数内部生成。从 SAR 切换到 SDM 后模型待估参数从 k 个变成 2k 个自由度消耗增加。样本量不够大的时候SDM 的估计方差会显著变大。如果 T 很小比如只有 5 年我建议先跑 SAR 和 SEM 作为参照再决定是否需要 SDM 的复杂度。是否选择 SDM可以用后面第 4 章讲的 LR 检验去看 WX 联合显著性。4. 输出结果怎么读结构体字段、直接与间接效应、模型选择检验4.1 结果结构体里到底藏了哪些字段Elhorst 代码的输出统一封装在结构体 results 里字段名有规律。搞清楚每个字段的含义比反复翻代码更高效。下表是固定效应模型最常见的字段说明。字段含义使用场景beta解释变量系数向量核心估计结果tstat各系数的 t 统计量显著性判断rho空间自回归系数空间依赖强度lambda空间误差系数SEM 用误差项空间相关sige误差方差估计模型拟合质量lik对数似然值模型比较yc / xcdemean 后的 y 和 x边际效应计算direct / indirect / totalSDM 的效应分解解释变量边际影响这里我想提醒一点beta 的 t 统计量是基于渐近分布计算的面板数据 T 较小或者 N 较小时t 值会偏乐观。论文报告时最好同时给出 bootstrap 标准误做稳健性核查不要只依赖 tstat 字段。4.2 直接效应与间接效应什么情况下必须看如果一个模型里没有 WX 项比如纯 SAR解释变量的边际效应就是 beta 本身直接读系数即可。但 SDM 包含 WX 项之后某个解释变量 x_k 的边际影响会分两条路径传导本地区 x_k 变化对本地 y 的影响直接效应以及本地 x_k 变化通过 WX 传导到邻居、邻居再反馈回本地的闭环影响间接效应。SDM 里 WX 的系数 theta 不能直接解读为“邻居影响”必须看分解后的效应。% 从SDM结果中提取效应做报表 direct_effect results_sdm.direct; % k×1 向量 indirect_effect results_sdm.indirect; % k×1 向量 total_effect results_sdm.total; % k×1 向量 % 效应对应的解释变量名称 var_names {x1, x2}; for i 1:length(var_names) fprintf(%s: 直接效应%.4f, 间接效应%.4f, 总效应%.4f\n, ... var_names{i}, direct_effect(i), indirect_effect(i), total_effect(i)); end效应分解的理论基础是偏导数矩阵。在 SDM 模型里对 y 求关于 x_k 的偏导会得到一个 N×N 的矩阵主对角线元素取平均就是直接效应非对角线元素按行或按列取平均就是间接效应。Elhorst 代码返回的 direct、indirect、total 已经完成了这些平均你直接拿去写表即可。有一个很容易被忽略的细节间接效应可能和直接效应符号相反这在实证里不少见。比如本地增加某政策变量直接效应为正但邻居受到负面溢出整体总效应存疑。写论文时遇到这种情况一定要分路径讨论不能只报总效应了事。4.3 LM、LR、Wald 三个检验怎么配合使用模型选择是空间面板实证绕不开的环节。Elhorst 包提供了三类检验用途各有侧重配合使用能覆盖从模型设定到稳健性的完整链条。LM 检验用于判断空间依赖的存在形式分为空间滞后检验lmlag_panel和空间误差检验lmerror_panel。常见做法是先跑一个不含空间项的普通面板回归固定效应 OLS然后对残差做 LM 检验。如果两种 LM 统计量都显著大概率需要 SDM只有滞后显著选 SAR只有误差显著选 SEM。% 先跑普通固定效应模型作为基准 % 注意: Elhorst包内没有单纯的面板OLS函数, 常见做法是手动demean后用OLS % 这里以包内已有的检验函数为例 [LM_lag, pval_lag] lmlag_panel(y, x, W, T); fprintf(LM空间滞后检验: LM%.4f, p%.4f\n, LM_lag, pval_lag); [LM_err, pval_err] lmerror_panel(y, x, W, T); fprintf(LM空间误差检验: LM%.4f, p%.4f\n, LM_err, pval_err);LR 检验通常用于嵌套模型的比较比如判断 SDM 是否可以退化为 SAR。SDM 退化为 SAR 的约束条件是所有 theta 系数联合为零lratios 函数可以给出这个约束的 LR 统计量。Wald 检验的用途类似但 LR 需要估计约束模型和无约束模型两个结果Wald 只需要无约束模型计算更省。我的习惯是四步走先 LM 检验判断空间项形式再跑 SAR 和 SDM 两个模型用 LR 检验看 WX 联合显著性最后对最终选定的模型做 Wald 稳健性核查。这个流程下来审稿人常见的“模型设定依据不足”质疑基本能堵住。5. 避坑记录Elhorst 代码包最常见的五类翻车现场5.1 权重矩阵维度写错报错却指向数据现象调用 sar_panel_FE 时报错 “Matrix dimensions must agree”或者程序直接卡死内存占用飙升。这时候你往往会去检查 y 和 x 的维度但它们没问题。原因权重矩阵 W 应该是 N×N代码内部会自动把它扩展到 NT×NT 参与运算。如果你误以为 W 需要传 NT×NT 的完整块对角矩阵等于把矩阵扩展后的维度又翻了几倍内存直接爆掉不常见的还会出现维度错位报错。解决确认 W 是 N×N 而不是 NT×NT。用 size(W) 看一眼再把 W 转成 sparse 稀疏矩阵存储。Elhorst 函数内部接受稠密矩阵但高维时稀疏矩阵能省一半内存。5.2 面板排序错位结果看似正常实则全错现象模型正常收敛rho 落在 0.3 到 0.6 之间t 值也显著但 beta 的符号和理论预期相反或者数值大得离谱。原因y 和 x 的个体排列顺序不一致。常见于数据清洗时删除了某个个体的部分时期导致后续个体的时间序列整体前移而你没有重新排序。解决进模型之前用上一章的 sortrows 方法做强制排序再抽两个个体核对行列对应关系。从那以后我每次处理新面板数据都会先做一次排列核对再进估计。5.3 rho 贴着边界跑多半是 W 标准化出了问题现象rho 跑到 0.98 甚至 0.999似然值却还在上升迭代到边界才停下。这种情况下结果极不稳定去掉几个样本 rho 就大变。原因W 行标准化之后谱半径接近 1空间滞后项的系数理论上限也随之贴近 1。如果 W 本身构造有误比如存在孤岛、行和过大rho 会被顶到边界。另一种可能是 W 里包含了对角线未清零的冗余连接。解决检查 max(abs(eig(W))) 是否约等于 1检查 W 的对角线是否为全零检查是否有某一行全为 0孤岛个体。修正后重新估计rho 通常会回到合理区间。5.4 rar 解压失败与 Matlab 路径里的中文现象双击 rar 包提示解压错误或文件损坏强行解压后部分 .m 文件缺失。还有人把文件放在带中文的路径下Matlab 能找到文件但运行时报错“未定义函数或变量”。原因rar 下载不完整、压缩工具兼容性问题或者压缩包本身有密码保护。中文路径导致 Matlab 的函数解析失败是另一个独立问题跟代码无关但很容易被误判为代码有 bug。解决优先用 7-Zip 打开 rar选择“提取到当前文件夹”。如果提示文件头损坏换 WinRAR 的修复功能工具 → 修复压缩文件再提取。解压后把整个文件夹拷到纯英文路径下比如 D:\elhorst_model_new再添加进 Matlab 路径。5.5 固定效应和随机效应跑出来几乎一样现象sar_panel_FE 和 sar_panel_RE 给出的 beta 相差不到 0.01rho 也基本一致看似模型稳健实际可能是数据里的个体异质性太弱或排序被破坏。原因固定效应模型通过 demean 消除了个体均值随机效应模型假设个体效应与解释变量不相关。如果真实数据里个体效应本身不显著FE 和 RE 的估计自然趋于一致。另一种可能是数据在构造时把时间均值也一并去掉了相当于人为削弱了个体差异。解决先看 RE 模型输出的 etarho个体效应占比如果接近 0说明个体效应确实弱结果一致是正常的不用纠结。如果是数据先做了全局去均值导致的回到原始数据重新处理。我一般会在估计前跑一个 F 检验判断个体效应是否显著再决定是否值得对比 FE 和 RE。6. 进阶技巧用模拟数据验证这套代码的估计精度拿到新版本的 Elhorst 代码包我不会直接拿来跑真实数据而是先走一遍蒙特卡洛模拟。空间计量估计器很容易在小样本下产生偏差用已知参数生成的数据检验代码的正确性比在真实数据上反复试错高效得多。% 蒙特卡洛模拟: 验证 sar_panel_FE 能否恢复真实参数 rng(42); % 固定随机种子, 保证结果可复现 N 100; % 个体数 T 10; % 期数 k 2; % 解释变量个数 rho_true 0.5; % 真实空间自回归系数 beta_true [1; -0.8]; % 真实解释变量系数 % 构造行标准化的 W W zeros(N, N); for i 1:N % 每个个体和最近的5个邻居相连, 模拟空间邻接 for j 1:5 nb mod(i j, N); if nb 0, nb N; end W(i, nb) 1; end end W normw(W); % 生成面板数据 x1 randn(N*T, 1); x2 randn(N*T, 1); x [x1, x2]; % 向量化生成空间面板数据 % y (I - rho*W)^(-1) * (X*beta epsilon) I_NT kron(eye(T), inv(eye(N) - rho_true * W)); eps randn(N*T, 1) * 0.5; y I_NT * (x * beta_true) eps; % 用包内函数估计 results sar_panel_FE(y, x, W, T); % 对比真实值与估计值 fprintf(rho: 真实%.2f, 估计%.4f\n, rho_true, results.rho); fprintf(beta1: 真实%.2f, 估计%.4f (偏差%.4f)\n, ... beta_true(1), results.beta(1), results.beta(1) - beta_true(1)); fprintf(beta2: 真实%.2f, 估计%.4f (偏差%.4f)\n, ... beta_true(2), results.beta(2), results.beta(2) - beta_true(2));模拟数据的构造逻辑是反解空间面板模型先设定真实的 rho 和 beta利用矩阵求逆生成 y再用代码包估计对比估计值和真实值之间的偏离。偏差在 0.1 以内说明代码实现可靠偏差过大则需要回到 W 矩阵或数据排序上排查。参数选择上有几个值得注意的地方。N 和 T 的设定直接关系偏差程度T 固定为 10 时N 从 50 提高到 200rho 的估计偏差会明显缩小这正是空间计量里的常见规律。epsilon 的方差设为 0.5是为了模拟中等信噪比的数据避免噪声过小导致估计看起来过于完美。我用这套方法踩过不止一次坑。某个网传的代码版本在 SDM 模型的效应分解上存在符号反转 bug就是靠模拟数据比对真实值抓出来的。从那以后我每次拿到新版本代码都强制自己先跑一遍模拟再上真实数据这个习惯帮我省掉了大量实证阶段反复排查的时间。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
