简介面向空间计量经济学研究者与高年级学生的动态空间杜宾模型MATLAB实现包专注空间动态回归建模支持时空滞后效应分析区别于传统静态模型可用于研究区域经济增长、产业集聚等面板数据中的空间溢出与时间动态特征。压缩包共19个文件以12个M脚本为主涵盖模型估计主函数含SAR、SARAR变体、配套示例与绘图输出脚本另含4个xlsx数据文件及1份PDF说明文档整体仅337KB便于下载与快速上手。已有1948人学习使用资源内含可直接运行的示例脚本、空间权重矩阵构造工具及生产率与集聚相关数据读者对照PDF笔记逐步运行可掌握动态空间面板模型的估计流程、结果解读与实战应用。适合需要开展空间计量实证研究或扩展论文模型的中高级用户。1. 动态面板空间杜宾模型空间计量里最值得下功夫的一类模型动态面板空间杜宾模型Dynamic Panel Spatial Durbin Model简称动态 SDM是静态空间杜宾模型在时间维度上的扩展它同时在方程右侧加入被解释变量的一阶时间滞后、空间滞后以及时空滞后项使得模型能同时捕捉时间惯性、空间溢出效应以及二者的交互影响。这个资源包给的不是空壳代码而是一整套能在 MATLAB 里直接跑通的动态面板 SDM 实现主程序 sar_jihai_time.m、似然函数 f2_sar_jihai_time.m、空间权重矩阵生成工具 wrook.m外加一份说明动态面板极大似然估计原理的笔记 note-sdpd-mle.pdf。配套数据是某地区高技术产业生产率与集聚面板数据非常适合正在做空间计量实证、写学位论文或投稿期刊的从业者也适合刚入门空间面板、想对比静态模型与动态模型差异的初学者直接上手复现。2. 为什么必须用动态空间杜宾模型静态模型在时间维度上的三个先天短板2.1 静态 SDM 只能给出当期均衡关系解释不了惯性静态空间杜宾模型的经典设定是y ρWy Xβ WXθ ε其中 ρ 是空间自回归系数Wy 是空间滞后项WX 是解释变量的空间滞后θ 是空间溢出系数。静态模型本质上假设所有变量在同一期内同时决定回归结果描述的是均衡状态下的空间交互关系。但现实中的经济数据几乎都有惯性——上一期的产出水平会深刻影响本期产出上期的技术溢出也需要时间才能在空间上扩散。静态模型把时间项塞进扰动项里导致两个直接后果一是遗漏解释变量造成的内生性偏差二是如果真实数据生成过程包含时间滞后静态模型的系数估计就是非一致的。动态 SDM 在方程右侧显式加入 yₜ₋₁把时间惯性从扰动项里捞出来。以资源包中主程序 sar_jihai_time.m 对应的模型形式为例yₜ τyₜ₋₁ ρWyₜ ηWyₜ₋₁ Xₜβ WXₜθ μ εₜ其中 τ 是时间滞后系数衡量惯性ρ 是当期空间溢出η 是时空滞后系数衡量的是上一期邻居的 y 对本期本地区 y 的影响。这个设定比静态模型多出了三个维度的信息时间惯性、空间溢出、以及两者交叉形成的时空扩散路径。2.2 静态模型估计的空间溢出效应被平均化动态模型能区分短期与长期效应静态 SDM 的空间溢出系数 θ 是一个打包值它无法区分溢出是当期发生的还是跨期累积的。例如地区 A 的研发投入提高 1%对地区 B 产出的影响可能第一年为 0.3%第二年积累到 0.8%静态模型只能估计出一个两者平均后的数。动态面板 SDM 则能从估计出的 τ、ρ、η 中解析出短期直接效应、短期间接效应、长期直接效应、长期间接效应四个量。这对于政策评价极其关键——短期效应回答当年见效多少长期效应回答稳态下累计见效多少。资源包中 note-sdpd-mle.pdf 的推导部分核心就是在讲动态条件下极大似然函数中雅可比项的处理。因为静态 SDM 只需处理一次 ln|Iₙ − ρW|而动态 SDM 每个时间截面都要处理 ln|Iₙ − ρW|同时还要考虑初始值 y₁ 的分布设定。这个问题处理不好估计结果就是有偏的后面的效应分解全都建立在错误估计量上。2.3 普通 OLS 和静态 ML 在动态空间设定下完全失效空间面板模型中Wyₜ 是内生的因为它与 εₜ 相关yₜ₋₁ 在加入个体固定效应 μ 后也是内生的Nickell 偏误而这种偏误在空间依赖下会被 ρ、η 进一步放大。OLS 估计动态面板空间模型时τ 的偏误不会随着 N 增大而消失只会在 T 增大时减弱。这就是为什么资源包宁可写极大似然——虽然 MLE 对分布假设敏感但在小样本和典型空间面板设定下它仍然比 GMM 更容易收敛、更稳定。sar_jihai.m 是静态版本sar_jihai_time.m 是动态版本两个程序并排放在包里最适合做对比把同一份数据分别跑静态和动态看系数和显著性发生了哪些变化这是理解动态模型价值的最快捷路径。3. 跑通第一个动态面板空间杜宾模型从压缩包到回归结果的全流程3.1 压缩包内文件结构与功能定位拿到压缩包后先别急着运行花十分钟把文件分好类。这里给出我拆包后的文件对应关系文件/目录功能角色说明sar_jihai.m静态 SAR 主程序不含解释变量 X 的纯空间自回归模型适合做基准对照sar_jihai_time.m动态 SAR 主程序含 yₜ₋₁、Wyₜ、Wyₜ₋₁ 的动态空间自回归模型f_sar_jihai.m静态 SAR 似然函数被 sar_jihai.m 调用的目标函数返回负对数似然值f2_sar_jihai_time.m动态 SAR 似然函数动态模型的 MLE 目标函数关键改动集中在这个文件SARAR.mSARAR 模型主程序空间自回归带空间自相关误差项的模型用于误差项存在空间依赖时的备选估计prt_sardynamic.m动态结果输出工具格式化打印动态模型的系数、t 统计量、拟合优度prt_sp.m / prt_reg.m通用结果输出工具宏观经济计量工具包中的通用打印函数wrook.m空间权重矩阵生成器根据经纬度或邻接关系生成 Rook 一阶邻接权重矩阵note-sdpd-mle.pdf方法笔记动态面板数据 MLE 推导笔记建议在跑代码前通读一遍data.xlsx面板数据被解释变量与解释变量的面板数据按地区-年份排列matrix.xlsx空间权重矩阵数据地区间空间关系的原始数据需要读入后用 wrook.m 处理实际运行中main 入口建议使用 sar_jihai_example.m因为它已经把数据读取、权重矩阵生成、模型估计、结果打印的完整流程串好了可以避免自己拼 matlab 脚本时搞乱数据顺序。3.2 示例脚本跑通完整估计以 sar_jihai_example.m 为入口在 MATLAB 中把当前路径切换到解压目录后最直接的复现方式是写一个驱动脚本按照以下逻辑执行。这里给出一个标准模板%% 清理工作区并设置路径 clear; clc; addpath(你的解压路径); %% 读取面板数据与空间矩阵 data readmatrix(data.xlsx); % 第一列是地区ID第二列是年份其余列是变量 Wraw readmatrix(matrix.xlsx); % 原始空间邻接矩阵非标准化形态 [n, T] size(data); % n为地区数T为年份数 y reshape(data(:, 3), n, T); % 假设第3列是被解释变量重塑成 n×T X reshape(data(:, 4:end), n, T, size(data, 2) - 3);这里有一个重要设计面板数据在 Excel 中是长格式每个地区-年份组合一行但 MATLAB 估计程序需要的是宽格式每个变量是 n×T 的矩阵。reshape 操作正是完成这一步的。如果你的数据顺序稍有不同务必先检查 data.xlsx 的列顺序否则后续结果全部错位。生成空间权重矩阵%% 生成标准化的 Rook 一阶邻接权重矩阵 W wrook(Wraw); % 输入原始0-1邻接矩阵输出行标准化后的W Wsp sparse(W); % 转稀疏矩阵提高后续矩阵运算效率wrook.m 的逻辑是先判断两个地区是否拥有公共边界是则赋 1否则赋 0随后对每一行做行标准化使每行元素之和为 1。行标准化后的权重矩阵可以保证空间滞后项 Wy 的经济含义为邻居的加权平均。接下来调用动态模型估计%% 动态SAR模型估计 info.lflag 0; % 0精确对数行列式计算1近似计算 info.model 1; % 1个体固定效应2随机效应3时间固定效应 info.lndet []; % 若已有预计算的lndet可传入加速 result sar_jihai_time(y, X, W, info); %% 打印结果 prt_sardynamic(result);参数说明info.lflag 控制对数行列式 ln|Iₙ − ρW| 的计算方式。精确计算在 n 400 时完全可行当 n 超过 500 或矩阵频繁迭代时建议改为 info.lflag 1 走蒙特卡洛近似速度能提升数倍但精度会轻微下降。info.model 决定固定效应形式个体固定效应是空间面板中使用最广泛的选择。prt_sardynamic 会输出 τ、ρ、η 的估计值及对应 t 统计量还会报告 log-likelihood、R² 等信息。%% 对照同样数据跑静态SAR模型 info_static info; result_static sar_jihai(y, X, W, info_static); prt_sp(result_static);对比两次输出的 Log-likelihood 和系数你能直观看到引入动态项后模型拟合是否显著改善。这是论文中最常用来论证动态模型优于静态模型的关键证据。3.3 权重矩阵的生成细节直接影响所有估计结果的下游参数权重矩阵是空间计量里最玄学的部分——同样的数据换一种权重矩阵ρ 和 η 的显著性可能完全翻转。wrook.m 生成的是一阶 Rook 邻接矩阵适合基于地理邻接关系的产业集聚问题。如果你处理的是经济距离或引力模型还需要改代码。这里给出 wrook.m 的核心逻辑function W wrook(Wraw) % Wraw: 0-1邻接矩阵Wraw(i,j)1表示地区i与地区j相邻 n size(Wraw, 1); W Wraw; % 行标准化 for i 1:n s sum(W(i, :)); if s 0 W(i, :) W(i, :) / s; end end end这个函数逻辑很简单但有几个容易踩的细节第一Wraw 对角线必须是 0自己不能算自己的邻居第二孤立地区某一行全为 0不会被标准化对应行保持全 0表示它没有邻居。如果你的数据样本里包含某个无邻接的地区直接跑动态模型可能会出现 ρ 估计值异常因为该地区的空间滞后项恒为 0等同于一个异常值注入。建议在跑模型前先检查 sum(Wraw, 2) 是否出现 0 值。4. 核心程序拆解sar_jihai_time.m 与 f2_sar_jihai_time.m 的估计逻辑4.1 动态空间自回归模型的 MLE 目标函数怎么写理解动态 SDM 的估计逻辑关键在于看 f2_sar_jihai_time.m 这个目标函数。动态模型的对数似然函数在经典静态基础上增加了时间滞后项的雅可比调整整体形式为lnL −(nT/2)ln(2πσ²) T·ln|Iₙ − ρW| − (1/2σ²)·Σₜ eₜ eₜ其中 eₜ yₜ − τyₜ₋₁ − ρWyₜ − ηWyₜ₋₁ − Xₜβ − WXₜθ − μ。代码实现时为了减少矩阵运算量会先用 Frisch-Waugh-Lovell 定理把固定效应 μ 消去再对转换后的数据做迭代搜索。sar_jihai_time.m 的主体是两层循环外层循环扫描候选 ρ 值通常是在 (−1, 1) 区间内做 100 个格的网格搜索内层用最小二乘快速计算给定 ρ 下的 β、τ、η、θ 和 σ²。最后取使似然值最大的 ρ 作为最优估计。function llike f2_sar_jihai_time(param, y, x, W, T, n, info) % param: [rho; beta; tau; eta; sigma2] 的列向量 rho param(1); beta param(2:1size(x,2)); tau param(end-3); eta param(end-2); sigma2 param(end); In speye(n); A In - rho * W; % 空间滤波矩阵 yt y(:); % 将 y 拉成列向量 % 构建 H y_t - tau*y_{t-1} - rho*W*y_t - eta*W*y_{t-1} % 这里需要按时间切片循环处理 e zeros(n*T, 1); for t 2:T y_prev y(:, t-1); y_cur y(:, t); Wy_cur W * y_cur; Wy_prev W * y_prev; X_cur x(:, :, t); e((t-1)*n1 : t*n) A * y_cur - tau * y_prev - eta * Wy_prev ... - X_cur * beta; end % 对数似然值 llike -(n*(T-1)/2)*log(2*pi*sigma2) (T-1)*log(det(A)) ... - (1/(2*sigma2)) * (e * e); llike -llike; % 返回负值便于最小化 end代码逻辑说明目标函数接收参数向量 param取出 ρ、β、τ、η、σ² 后构造残差向量 e。核心是循环中 A*y_cur − τ*y_prev − η*Wy_prev 这三项分别对应当期空间滤波、时间惯性、时空滞后。注意这里没有直接写 WXθ 项说明动态 SAR 模型和动态 SDM 在代码实现上有区别——SDM 还需要构造 WX 并追加到解释变量矩阵中sar_jihai_time.m 名称里的 sar 表明它当前实现的是纯空间自回归结构如果你需要用动态 SDM需要手动把解释变量的空间滞后列拼接到 X 中。4.2 初始值与迭代收敛极大似然搜索中最容易被低估的环节MLE 迭代的本质是在似然曲面上找峰但这个曲面在 ρ 接近边界时非常平坦导致搜索算法可能在 ρ 0.95 和 ρ 0.98 之间反复震荡。资源包中的 sar_jihai_time.m 采用了两阶段策略第一阶段在 [−0.99, 0.99] 区间均匀取 30 个初始点对每个点做一次完整迭代选择似然值最高的结果作为最终迭代的初值第二阶段使用 fminsearch 或 fminbnd 做精细搜索。这段逻辑实现如下% sar_jihai_time.m 内部的两阶段搜索示意 rmin -0.99; rmax 0.99; ngrid 30; rgrid linspace(rmin, rmax, ngrid); llmax -inf; for i 1:ngrid rho_init rgrid(i); % 用当前rho初值计算条件最大似然 [ll, ~] f2_sar_jihai_time([rho_init; beta0; tau0; eta0; sigma2_0], ... y, x, W, T, n, info); if ll llmax llmax ll; rho_best rho_init; end end % 用rho_best作为最终搜索起点 param0 [rho_best; beta0; tau0; eta0; sigma2_0]; options optimset(Display, iter); param_hat fminsearch((p) f2_sar_jihai_time(p, y, x, W, T, n, info), ... param0, options);参数说明ngrid 是网格搜索密度n 超过 200 时建议降到 20否则每次网格点都调用一次似然函数计算成本会很高。fminsearch 是 MATLAB 自带的 Nelder-Mead 算法对不可导或含平坦区域的函数比较鲁棒。如果你发现迭代半天都不收敛可以先检查目标函数是否返回了 NaN——这通常发生在线性代数运算中的矩阵接近奇异时也就是 ρ 接近 1/W 最大特征值的倒数。4.3 从静态到动态的改动对照三行核心代码的区别把 sar_jihai.m 与 sar_jihai_time.m 并排对比差异非常集中。静态模型的似然函数中残差构造为e (Iₙ − ρW)y − Xβ动态模型则多出两步一是滞后项 yₜ₋₁ 进入解释变量集合二是空间滞后项 Wyₜ₋₁ 也要进入。这意味着如果你仅仅在 sar_jihai.m 的代码里加一列滞后 y是远远不够的——还必须同步加入 Wyₜ₋₁并在空间滤波矩阵 A 的处理上区分当期项与滞后项。很多初学者只加 yₜ₋₁然后发现 η 估不出来或者显著性极差原因就在这里。正确的扩展方式是构造一个新的解释变量矩阵 Z [yₜ₋₁, Wyₜ₋₁, Xₜ, WXₜ]然后估计yₜ ρWyₜ Zδ μ εₜ。sar_jihai_time.m 内部实际做的就是这件事f2_sar_jihai_time.m 的残差构造循环体现了这一点。5. 避坑指南动态面板空间杜宾模型最常见的五个翻车点5.1 现象τ 估计值接近 1模型不收敛或结果异常原因被解释变量存在强烈的单位根或近单位根过程。动态面板模型要求 |τ| 1 以保证稳定性如果变量本身是 I(1) 过程例如未经处理的产出水平值τ 会无限接近 1模型在边界处无法收敛。解决先对被解释变量做单位根检验必要时取对数差分或增长率形式。产业集聚类数据通常用区位熵或密度指标这些指标本身是平稳的但如果直接用原始产值极大可能翻车。另一个处理办法是加入时间趋势项或者改用 longterm.m 对应的长期模型设定来缓解。5.2 现象ρ 总是被推到 0.99 以上且网格搜索图完全不呈现单峰原因权重矩阵 W 与经济距离不匹配。Rook 邻接矩阵只考虑地理边界相接如果研究的是产业间技术溢出地理邻接可能根本无法有效刻画溢出渠道。此时 ρ 会被推高以弥补权重矩阵解释力度不足的问题。解决更换权重矩阵类型。常见的做法是构造经济距离权重矩阵——用地区间人均 GDP 差异的倒数作为权重或者构造引力模型权重——用 GDP 乘积除以地理距离。更稳妥的做法是同时跑三套权重Rook、Queen、经济距离做敏感性分析在论文中报告结果是否稳健。5.3 现象η 时空滞后项不显著但理论上明显应该存在原因面板数据的年份间隔过大。如果 T 5 且每期间隔 5 年时空滞后效应可能已经在期內完全衰减统计上捕捉不到。另一个常见原因是数据排列顺序错误——MATLAB 中如果 yₜ₋₁ 取成了同一年的另一列数据η 自然不显著。解决检查数据排列。确认 data.xlsx 中每个地区的年份是严格连续的并且 reshape 之后 y(:, t) 对应的是第 t 年。建议在估计前画出 y 的热力图肉眼检查空间分布模式是否随时间推移出现明显的波瓣状扩散如果有说明时空滞后确实存在η 不显著的问题在数据质量或权重矩阵上。5.4 现象静态模型和动态模型的系数符号完全相反原因动态模型中 τyₜ₋₁ 吸收了相当一部分原来静态模型中 X 的解释力。如果解释变量本身具有较强的惯性例如研发投入历年变化不大静态模型中 X 的系数包含了惯性成分动态模型把惯性剥离给 τ 后X 的系数可能缩小甚至变号。这是正常现象不是代码 bug。解决在论文里解释清楚——静态模型估计的是总效应动态模型估计的是净效应。惯例是报告两套结果并重点解释动态模型净效应的经济含义。如果动态模型中 β 由显著变不显著说明该解释变量的影响主要通过惯性传导而非当期直接作用这本身就是一个有价值的结论。5.5 现象运行时报错 Insufficient number of observations 或维度不匹配原因动态模型需要滞后项 yₜ₋₁所以实际使用的观测是 T−1 个时间截面而不是 T 个。如果程序中 y 和 X 的维度没有相应调整就会出现矩阵乘法维度对不上。解决检查输入到 f2_sar_jihai_time.m 的 y 是否已经提前丢弃了第一年数据。有经验的写法是y y(:, 2:end); % 丢弃第一期因为需要 y_{t-1} 构造滞后 X X(:, :, 2:end); T T - 1;然后重新估计。这一步遗忘是初学者报错的第一大原因静态模型没有这个要求所以从静态代码改动态时特别容易忽略。6. 结果验证与进阶用法动态 SDM 的效应分解与稳健性自检6.1 先算直接效应与间接效应再谈系数解释动态空间杜宾模型的系数 β 不能直接解释为边际效应因为空间溢出通过反馈循环会回到本地区。正确的做法是基于估计出的 τ、ρ、η、θ 求取偏导数矩阵。对于动态 SDM第 t 期的短期直接效应是矩阵 (Iₙ − ρW)⁻¹(β Wθ) 的对角线均值短期间接效应是对角线外元素的行均值。长期效应则要在短期基础上除以 (1 − τ) 的调整因子因为长期中时间惯性会持续放大影响。在 MATLAB 中实现短期与长期效应分解%% 基于估计结果的效应分解 rho result.rho; beta result.beta; theta result.theta; tau result.tau; A_inv inv(eye(n) - rho * W); % 短期直接效应与间接效应 short_mat A_inv * (beta theta); % 这里假设X是单变量多变量时需逐变量构造 direct_short mean(diag(short_mat)); indirect_short mean(sum(short_mat, 2) - diag(short_mat)); % 长期效应除以 1-tau 调整 long_mat short_mat / (1 - tau); direct_long mean(diag(long_mat)); indirect_long mean(sum(long_mat, 2) - diag(long_mat)); fprintf(短期直接效应: %.4f\n, direct_short); fprintf(短期间接效应: %.4f\n, indirect_short); fprintf(长期直接效应: %.4f\n, direct_long); fprintf(长期间接效应: %.4f\n, indirect_long);这段代码说明效应分解的输出比系数本身更有说服力。如果间接效应显著而直接效应不显著说明该变量主要通过空间溢出路径影响邻地。论文中报告这四个数字比单独罗列 β 更有冲击力。6.2 稳健性检验的三个必做动作做稳健性检验时我一般强制自己走三关。第一关是换权重矩阵——把 Rook 换成 Queen 邻接或者经济距离权重观察 ρ、η 显著性是否保持。第二关是换估计方法——用 SARAR.m 估计带空间自相关误差项的模型看核心结论是否发生变化。第三关是改变动态项设定——尝试只含 yₜ₋₁ 不含 Wyₜ₋₁ 的模型与完整动态模型做似然比检验。三个动作做完还稳健的结论审稿人基本不会再在这上面做文章。6.3 最终建议跑模型前先跑一遍数据诊断我现在每次拿到一份新的空间面板数据都强制自己先花半小时做三件事画被解释变量的时间趋势图看是否存在明显的共同趋势算一下变量的组内相关系数判断固定效应与随机效应的选择对空间权重矩阵做特征值检查确保 ρ 的有效参数空间是 (−1/λ_min, 1/λ_max)。这三步看起来朴素但能挡掉一大半后面的不收敛和假显著问题。这些年拆过的空间计量代码包里这套动态 SDM 的 MATLAB 实现是我见过少有的能直接跑通、又不靠黑箱命令的工具。它的价值不仅在代码本身更在于那份 note-sdpd-mle.pdf 把动态面板模型的似然推导写得足够清楚照着推一遍再回来看程序每个函数都能对得上号。希望你也能在这套代码上跑出自己的结果希望这次拆解帮你少走几个月的弯路。本文还有配套的精品资源点击获取