简介动态面板空间杜宾模型资源包聚焦空间计量经济学中的动态面板回归面向研究空间溢出效应的经济学、地理学及区域科学领域人士尤其适用于需要从传统静态模型转向动态框架的实证研究者。该资源共包含19个文件以12个Matlab脚本为核心覆盖模型估计、结果输出与动态效应计算另配4个Excel数据文件和1份PDF理论说明压缩包整体约337KB。目前已有1949人浏览学习适合硕博论文、课程作业及期刊复现等场景。通过资源内置的SARAR、长期效应等程序读者可快速上手空间动态杜宾模型的设定与估计清晰对比其与传统静态模型的差异并结合示例数据直接运行验证。配套PDF文档则对极大似然估计方法进行补充阐释帮助深入理解动态空间面板建模的细节与陷阱。1. 动态面板空间杜宾模型当空间溢出碰上面板惯性静态 SDM 不够用了动态面板空间杜宾模型很多人简称动态 SDM是把“被解释变量的一期滞后”和“邻居当期的被解释变量”同时放进回归方程的面板空间计量模型。它解决的问题很具体在做区域政策评估或空间溢出研究时既要解释“本地上期状态对本期的影响”又要把“隔壁地区当期的互相反馈”分离出来此时只用静态空间杜宾模型会把时间路径依赖混进误差项得到偏高的空间滞后系数。这篇笔记写给要做实证、但不想在方法论上翻车的读者覆盖模型设定、权重矩阵、估计取舍、可复现的脚本和五类典型坑位末尾再给汇报论文时的效应分解技巧。2. 动态 SDM 的方程到底长什么样从静态到动态少加了哪一项2.1 完整方程与每一项的经济含义大多数实证研究者先跑熟的是静态空间杜宾模型写成y_it ρ * Σ_j w_ij * y_jt x_it * β Σ_j w_ij * x_jt * θ μ_i ε_it这个式子只处理了“当期邻居对我”的空间交互。把时间滞后项加进去后完整动态 SDM 写成y_it τ * y_i,t-1 ρ * Σ_j w_ij * y_jt x_it * β Σ_j w_ij * x_jt * θ μ_i γ_t ε_it其中 τ 是时间滞后系数衡量本地上期状态对本期的影响常被解释为惯性、路径依赖或黏性ρ 是空间滞后系数衡量邻居当期加权平均对本地的影响β 是解释变量自身的系数θ 是解释变量空间滞后项的系数反映邻居的 x 对本地 y 的溢出。μ_i 是个体固定效应γ_t 是时期效应。很多入门者容易踩的第一个坑是以为加了 y_i,t-1 就算“动态”了。严格来说动态 SDM 要求时间滞后、内生空间滞后Wy、外生空间滞后WX同时存在并且估计时必须处理 y_i,t-1 和 Wy_t 两个内生变量。如果只是把 y_i,t-1 扔进 xsmle 命令那不叫动态 SDM叫“静态 SDM 加了个回归量”。2.2 空间权重矩阵 W模型里最容易埋雷的部分空间权重矩阵 W 是 SDM 类和 SAR 类模型共用的核心输入。常见做法有三种邻接矩阵、距离衰减矩阵、K 近邻矩阵。选型上没有绝对标准但一定要在论文里说明选择理由。构造方式常见定义适用场景主要局限邻接矩阵共用边界或顶点取 1否则取 0行政区、地块、网格隔海相望或距离很远但有强联系时失效距离衰减w_ij 1 / d_ij^2 或 1 / d_ij城际溢出、经济地理距离幂次要人为设定结果对参数敏感K 近邻每个个体取距离最近的 K 个邻居样本区域密度不均匀时K 值选择有一定玄学需要做敏感性无论用哪种几乎都必须做行标准化w_ij* w_ij / Σ_j w_ij。行标准化之后W*y 这一项变成邻居 y 的加权平均ρ 的值被约束在可行区间里I - ρW也更容易保证可逆。两个最容易埋雷的细节一是对角线必须清零自己不能是自己的邻居二是距离矩阵如果不做行标准化谱半径可能超过 1后续求log|I - ρW|时直接得到 NaN。2.3 动态 SDM 与 SAR、SEM、静态 SDM 的选型关系模型的选取思路不是“越复杂越好”而是按照数据特征和理论问题一层层加上去。模型包含 y_t-1包含 Wy_t包含 WX典型估计策略回答什么问题SAR否是否ML邻居的 y 是否影响本地 ySEM否否否ML误差项是否存在空间相关静态 SDM否是是ML解释变量和被解释变量的双重溢出动态 SDM是是是QMLE / 空间 GMM短期效应与长期效应的区分我一般先跑静态 SDM用 LR 或 Wald 检验看能否退化为 SAR 或 SEM一旦发现被解释变量有明显的时间惯性比如上期污染物浓度对本期解释力很强再考虑动态 SDM。T 最好在 10 年以上T 太短时动态项和固定效应纠缠在一起估计结果极不稳定。3. 估计方法的取舍为什么 ML 在动态面板里会偏QMLE 与空间 GMM 怎么选3.1 内生性与 Nickell 偏误动态项的两个麻烦静态 SDM 常用的极大似然估计ML在动态面板里会遇到两个内生性问题。第一个来自 y_i,t-1个体固定效应 μ_i 进入所有期y_i,t-1 和 μ_i 必然相关这在经典动态面板里叫 Nickell 偏误。T 越小偏误越严重T 只有 10 到 20 期时τ 的估计会被明显向下拉。第二个问题来自 Wy_t空间滞后项是邻居当期的 y而邻居当期的 y 又受本地当期扰动影响Wy_t 和 ε_it 之间存在联立性也就是常说的“反射问题”。这两个问题叠在一起直接搬静态 ML 估计动态 SDM 的结果基本不能信。T 大到几十期时偏误才会逐渐消失但很多面板数据只有 10 到 15 年必须换估计策略。3.2 QMLE 偏误校正常见做法与算法骨架圈内最常见的做法是用带偏误校正的准极大似然估计QMLE。思路分三步先对面板做组内变换消掉个体固定效应再把空间滞后项和时间滞后项当作回归量用集中似然搜索 τ 和 ρ最后对组内变换导致的动态偏误做解析或自助法校正。很多现成的 Matlab 工具箱就是这么实现的。动手写自己的估计器之前至少先确认权重矩阵合格否则后面所有数值结果都是废的。我习惯在估计前跑一段矩阵检查% 检查行标准化后的空间权重矩阵 % W: N x N 空间权重矩阵要求对角线为零、行和为 1 maxEig max(abs(eig(full(W)))); if maxEig 0.9999 error(W 未正确行标准化谱半径接近或超过 1请检查对角线和对角线外的连接); end这段代码的逻辑是行标准化矩阵的最大特征值模理论上等于 1由于浮点计算允许一点点误差所以阈值设 0.9999。如果谱半径明显大于 1后面计算的log|I - ρW|在 ρ 较大时会直接变成无穷或 NaN。很多数据生成的流程跑着跑着出 NaN问题根本不在模型而在 W。3.3 什么时候换成空间系统 GMM如果 T 小于 15或者空间权重矩阵用的是经济距离这类与当期变量相关的矩阵QMLE 的偏误校正也不一定救得回来。此时可以考虑空间系统 GMM先做一阶差分消掉个体固定效应再用滞后两期及以上的水平项作为工具变量同时把 W*y 的空间滞后加入工具集。代价是工具变量数量容易爆炸。常见的翻车现场是 Sargan/Hansen 检验被拒绝原因往往不是模型错而是工具变量堆得太多或工具集里混入了当期内生项。空间系统 GMM 我一般只在 QMLE 结果明显异常、或者审稿人明确要求短面板稳健性时使用不会作为第一选择。4. 跑通动态面板空间杜宾模型的最小复现数据生成、权重矩阵与估计脚本4.1 准备一个可以自检的模拟数据DGP 先行不管读者是从别人手里拿到一个 .rar 压缩包还是攒了一套自己的面板数据第一步都是先让脚本能在一个“已知真实参数”的数据上跑通。Windows 下用 7-Zip 就能打开 rar 压缩包解压后先看里面有没有data、W、main三个部分如果没有就自己生成一份。我用一段可以完整运行的 Matlab/Octave 脚本生成动态 SDM 模拟数据。脚本里固定随机种子保证任何人复现的结果都一样rng(2024); % 固定随机种子保证结果可复现 N 60; T 10; % 60 个区域10 期面板 Wraw rand(N, N) 0.15; % 随机邻接关系平均每个区域约 9 个邻居 Wraw(1:N1:end) 0; % 对角线清零自己不是自己的邻居 W diag(1 ./ sum(Wraw, 2)) * Wraw; % 行标准化 Xmat randn(N, T); % 单个解释变量 mu randn(N, 1); % 个体固定效应 gam linspace(-0.5, 0.5, T); % 时期效应 tau 0.5; rho 0.3; % 真实参数 beta 1.0; theta 0.6; % 解释变量及其空间滞后系数 sig 0.3; % 扰动标准差 I speye(N); A I - rho * W; % 这一项必须可逆rho 要小于 1 Ymat zeros(N, T); for t 1:T xt Xmat(:, t); wxt W * xt; if t 1 rhs beta * xt theta * wxt mu gam(t) sig * randn(N, 1); else rhs tau * Ymat(:, t-1) beta * xt theta * wxt ... mu gam(t) sig * randn(N, 1); end Ymat(:, t) A \ rhs; end这段生成逻辑就是把动态 SDM 方程两边同时乘以(I - ρW)的逆把当期被解释变量显式解出来。t1 时没有上一期观测直接从静态分布启动之后每一期都同时受到上一期自己的影响、邻居当期的影响、解释变量和邻居解释变量的影响。参数说明tau0.5 表示时间惯性中等偏强rho0.3 表示空间溢出中等theta0.6 表示邻居的 x 对本地 y 有明显正向溢出。生成后建议先plot(Ymat(1:20,:))看一眼数据形态确认没有突变值。4.2 权重矩阵的行标准化与谱半径检查模拟数据里的 W 是随机生成的真实数据的 W 通常需要从地理坐标或行政区划构建。无论怎么构建估计前都按第 3 章的检查脚本跑一遍对角线是否为零、行和是否为 1、谱半径是否接近 1。如果发现某一行全为零说明存在孤立区域行标准化时会除零要提前处理成“以自身为邻居”或者删除该个体。% 检查是否有全零行 rowSum sum(Wraw, 2); if any(rowSum 0) warning(发现全零行请检查地理数据中是否有孤立区域); end全零行在邻接矩阵里很常见比如某个海岛县和任何陆地都不接壤。处理方式通常是两个要么用距离矩阵替代邻接矩阵要么把该区域从样本中剔除。强行保留会让行标准化那一步出 NaN后面的估计直接中断。4.3 QMLE 网格搜索脚本参数与初始值设置这里给出一版“集中似然 网格搜索”的完整脚本。它不是最前沿的偏误校正实现但把估计流程摆在了明面上对每个 τ 和 ρ先用外生变量回归得到 β、θ再计算残差和似然值最后取整个网格里似然最大的参数组合。真实研究中可以直接用现成的 QMLE 函数替换中间部分网格搜索的价值在于让人看清参数之间的耦合关系也方便排查为何优化器不收敛。% 输入Ymat(N x T)、Xmat(N x T)、W(N x N) % 输出tauHat、rhoHat、betaHat、thetaHat NT N * T; Yv Ymat(:); % 构造滞后项与空间滞后项 YlagMat [zeros(N, 1), Ymat(:, 1:end-1)]; WyMat W * Ymat; WXMat W * Xmat; % 组内去均值函数减去每个个体在 T 期上的均值 dm (z) z - kron(mean(reshape(z, N, T), 2), ones(T, 1)); Ydm dm(Yv); YlagDm dm(YlagMat(:)); WyDm dm(WyMat(:)); Xdm dm(Xmat(:)); WXDm dm(WXMat(:)); Z [Xdm, WXDm]; % 外生部分x 和 Wx % 参数网格先宽后密 tauGrid -0.5:0.05:0.9; rhoGrid -0.4:0.05:0.7; LL zeros(numel(tauGrid), numel(rhoGrid)); for a 1:numel(tauGrid) for b 1:numel(rhoGrid) tauVal tauGrid(a); rhoVal rhoGrid(b); % 构造剔除动态项后的因变量 ytmp Ydm - tauVal * YlagDm - rhoVal * WyDm; % 对剩余外生变量做 OLS得到 beta、theta bOLS Z \ ytmp; e ytmp - Z * bOLS; s2 (e * e) / (NT - size(Z, 2)); % 对数似然动态项贡献 空间雅可比项 logdet T * log(abs(det(I - rhoVal * W))); LL(a, b) -NT / 2 * (log(2 * pi * s2) 1) logdet; end end % 取网格最大值 [~, idxMax] max(LL(:)); [aIdx, bIdx] ind2sub(size(LL), idxMax); tauHat tauGrid(aIdx); rhoHat rhoGrid(bIdx); % 用最优 tau/rho 重新回归外生部分 ytmp Ydm - tauHat * YlagDm - rhoHat * WyDm; bOLS Z \ ytmp; betaHat bOLS(1); thetaHat bOLS(2); fprintf(tau_hat%.3f, rho_hat%.3f\nbeta_hat%.3f, theta_hat%.3f\n, ... tauHat, rhoHat, betaHat, thetaHat);代码里最关键的一行是ytmp Ydm - tauVal * YlagDm - rhoVal * WyDm它把时间滞后和空间滞后从因变量里“剔除”后剩余部分全部归因于外生变量。β 和 θ 的估计因此依赖于当前网格点的 τ 和 ρ这是集中似然的本质。网格范围要覆盖真实值的合理邻域动态面板 τ 一般落在 (-0.5, 0.9) 之间ρ 建议先参考静态 SDM 估计值再定。需要注意脚本保留了每个个体的第一期观测而严格研究中通常去掉第一期或者做初值处理否则滞后项的缺失值会影响 τ。4.4 从估计结果到短期/长期效应一张表说清楚空间杜宾模型里β 和 θ 不能直接当边际效应汇报。解释变量对 y 的影响会通过(I - ρW)^-1这个矩阵在网络里传播开所以要把直接效应、间接效应和总效应拆出来。单解释变量情况下短期效应矩阵是(I - ρW)^-1 * (βI θW)长期效应矩阵是((1-τ)I - ρW)^-1 * (βI θW)。% 从估计结果计算短期与长期效应 betaHat bOLS(1); thetaHat bOLS(2); Sshort (I - rhoHat * W) \ (betaHat * eye(N) thetaHat * W); Slong ((1 - tauHat) * I - rhoHat * W) \ (betaHat * eye(N) thetaHat * W); dirShort mean(diag(Sshort)); indShort mean(sum(Sshort, 2) - diag(Sshort)); dirLong mean(diag(Slong)); indLong mean(sum(Slong, 2) - diag(Slong)); fprintf(短期直接 %.3f间接 %.3f总效应 %.3f\n, ... dirShort, indShort, dirShort indShort); fprintf(长期直接 %.3f间接 %.3f总效应 %.3f\n, ... dirLong, indLong, dirLong indLong);这段代码的逻辑是短期效应矩阵只考虑当期空间反馈长期效应矩阵额外把时间累积算进去所以 τ 越大长期和短期的差距越大。汇报时可以整理成下面这种表审稿人最买账效应类型短期长期直接效应1.0231.214间接效应0.3870.562总效应1.4101.776注意不要把 βHat 填进“直接效应”一栏那是新手最容易犯的错。5. 动态 SDM 的五个常见坑现象、原因、解决一条条过5.1 似然值变成 NaN或者 ρ 被卡在网格边界现象网格搜索跑到一半log|I - ρW|输出 NaN或者最优 ρ 落在网格的最边缘比如 0.7 还想要更高。原因权重矩阵没有行标准化谱半径大于 1或者对角线没有清零自己和自己的空间滞后权重过强。解决先跑第 3.2 节的谱半径检查确认max(abs(eig(W)))在 1 附近。如果确实有个别行和不是 1重建矩阵不要手工缩放个别行直接统一行标准化。真实数据里如果存在全零行先处理孤立区域。5.2 滞后项串到别人的个体上τ 估计为负现象τ 的估计值和真实预期完全相反比如一个明显正相关的惯性过程估计出 -0.3并且 Wy 和 y_lag 的相关系数非常奇怪。原因面板数据从宽格式转长格式时个体排序没有按“个体-时间”排列滞后项[zeros(N,1), Y(:,1:T-1)]里的列顺序和当前行顺序错位前一期的数据跑到了别的个体上。解决进入估计前用sortrows或等价命令按(个体ID, 时间)排序排序后再构造YlagMat。我还会画前 20 个观测的连线图观察是否有断点跨过个体边界。5.3 把 β 系数直接当边际效应写进论文现象模型估计完正文说“解释变量每增加 1 单位被解释变量平均增加 β 个单位”。审稿人一眼看出问题。原因SDM 类模型通过(I - ρW)^-1把本地影响反馈到邻居再反馈回本地β 只是结构参数不是乘数结果。解决所有解释变量的汇报统一使用直接效应、间接效应和总效应不接受“先看系数再解释”。哪怕总效应方向和 β 相反也要以效应分解结果为准因为存在空间负反馈时 β 和总效应确实可能异号。5.4 动态项被当成“静态 SDM 加一个滞后变量”现象模型设定写了动态 SDM估计方法却填的是普通 ML且正文只讨论了 ρ 和 τ 的显著性没有讨论内生性处理。原因很多软件不能直接估计动态 SDM研究者就在静态模型里手动加y_t-1列用普通 ML 跑误以为加了滞后就是动态。解决务必在方法部分写明处理内生性的手段。如果是 QMLE要说明组内变换和偏误校正如果是空间系统 GMM要说明工具变量集。审稿人对“动态模型用静态估计”极其敏感这一条避不开。5.5 长期效应比短期效应大几个数量级现象长期间接效应是短期的几十倍甚至几百倍表格完全没法解释。原因(1 - τ)和ρ的关系出问题。长期效应矩阵((1-τ)I - ρW)接近奇异比如 τ0.85、ρ0.3 时最小特征值接近零逆矩阵被放大。解决估计后打印min(eig(full((1-tauHat)*I - rhoHat*W)))如果小于 1e-6说明模型接近非平稳。这时先看数据是否需要差分或者考虑把解释变量换成滞后项重新估计。不要强行解释一个条件数爆炸的矩阵。6. 汇报动态 SDM 结果的最后一步长期效应别只“乘个系数”要做矩阵稳健性6.1 从矩阵稳定性理解长期效应长期效应不等于短期效应乘以一个折扣。短期效应矩阵是(I - ρW)^-1 * (βI θW)它把当期空间传播算清楚了长期效应把时间累积也叠上去变成((1-τ)I - ρW)^-1 * (βI θW)。也就是说τ 通过改变I - ρW附近矩阵的条件数来影响长期结果而不是简单地把短期结果乘上1/(1-τ)。因此回归后的第一件事不是看系数表而是检查长期效应矩阵是否稳定。我习惯把下面这段判断放进所有动态 SDM 脚本里eigLong eig(full((1 - tauHat) * I - rhoHat * W)); if min(real(eigLong)) 1e-8 warning(长期效应矩阵接近奇异模型可能不平稳请检查 tau 和 rho 的组合); end这行命令比任何显著性检验都更能提前暴露问题。τ 和 ρ 单独看都显著合在一起却可能构成非平稳过程这类结果即使跑出来了也不能写进论文。6.2 权重矩阵敏感性换三种 W 看效应方向是否一致空间计量最容易被质疑的就是“换个权重矩阵结论就没了吗”。论文里一般要求做敏感性检验常见做法是准备三套矩阵邻接矩阵、距离倒数矩阵、K 近邻矩阵分别跑同一套估计脚本并输出三行效应。判断标准不是数值完全一致而是直接效应方向和总效应方向不要反转。如果邻接矩阵下溢出为正、距离矩阵下溢出为负说明结论对 W 过于敏感多半是模型设定漏了重要变量。我早期跑动态 SDM 时最喜欢盯着 τ 和 ρ 的显著性结果被审稿人一句“长期效应矩阵条件数太大”打回。后来养成一个习惯每次估计完先打印三个数——W 的谱半径、长期效应矩阵最小特征值、长期短期效应比值三个数都正常再开始写解释。希望这个习惯也能帮到你。本文还有配套的精品资源点击获取 动态面板空间杜宾模型实证指南：设定、估计与效应分解
简介动态面板空间杜宾模型资源包聚焦空间计量经济学中的动态面板回归面向研究空间溢出效应的经济学、地理学及区域科学领域人士尤其适用于需要从传统静态模型转向动态框架的实证研究者。该资源共包含19个文件以12个Matlab脚本为核心覆盖模型估计、结果输出与动态效应计算另配4个Excel数据文件和1份PDF理论说明压缩包整体约337KB。目前已有1949人浏览学习适合硕博论文、课程作业及期刊复现等场景。通过资源内置的SARAR、长期效应等程序读者可快速上手空间动态杜宾模型的设定与估计清晰对比其与传统静态模型的差异并结合示例数据直接运行验证。配套PDF文档则对极大似然估计方法进行补充阐释帮助深入理解动态空间面板建模的细节与陷阱。1. 动态面板空间杜宾模型当空间溢出碰上面板惯性静态 SDM 不够用了动态面板空间杜宾模型很多人简称动态 SDM是把“被解释变量的一期滞后”和“邻居当期的被解释变量”同时放进回归方程的面板空间计量模型。它解决的问题很具体在做区域政策评估或空间溢出研究时既要解释“本地上期状态对本期的影响”又要把“隔壁地区当期的互相反馈”分离出来此时只用静态空间杜宾模型会把时间路径依赖混进误差项得到偏高的空间滞后系数。这篇笔记写给要做实证、但不想在方法论上翻车的读者覆盖模型设定、权重矩阵、估计取舍、可复现的脚本和五类典型坑位末尾再给汇报论文时的效应分解技巧。2. 动态 SDM 的方程到底长什么样从静态到动态少加了哪一项2.1 完整方程与每一项的经济含义大多数实证研究者先跑熟的是静态空间杜宾模型写成y_it ρ * Σ_j w_ij * y_jt x_it * β Σ_j w_ij * x_jt * θ μ_i ε_it这个式子只处理了“当期邻居对我”的空间交互。把时间滞后项加进去后完整动态 SDM 写成y_it τ * y_i,t-1 ρ * Σ_j w_ij * y_jt x_it * β Σ_j w_ij * x_jt * θ μ_i γ_t ε_it其中 τ 是时间滞后系数衡量本地上期状态对本期的影响常被解释为惯性、路径依赖或黏性ρ 是空间滞后系数衡量邻居当期加权平均对本地的影响β 是解释变量自身的系数θ 是解释变量空间滞后项的系数反映邻居的 x 对本地 y 的溢出。μ_i 是个体固定效应γ_t 是时期效应。很多入门者容易踩的第一个坑是以为加了 y_i,t-1 就算“动态”了。严格来说动态 SDM 要求时间滞后、内生空间滞后Wy、外生空间滞后WX同时存在并且估计时必须处理 y_i,t-1 和 Wy_t 两个内生变量。如果只是把 y_i,t-1 扔进 xsmle 命令那不叫动态 SDM叫“静态 SDM 加了个回归量”。2.2 空间权重矩阵 W模型里最容易埋雷的部分空间权重矩阵 W 是 SDM 类和 SAR 类模型共用的核心输入。常见做法有三种邻接矩阵、距离衰减矩阵、K 近邻矩阵。选型上没有绝对标准但一定要在论文里说明选择理由。构造方式常见定义适用场景主要局限邻接矩阵共用边界或顶点取 1否则取 0行政区、地块、网格隔海相望或距离很远但有强联系时失效距离衰减w_ij 1 / d_ij^2 或 1 / d_ij城际溢出、经济地理距离幂次要人为设定结果对参数敏感K 近邻每个个体取距离最近的 K 个邻居样本区域密度不均匀时K 值选择有一定玄学需要做敏感性无论用哪种几乎都必须做行标准化w_ij* w_ij / Σ_j w_ij。行标准化之后W*y 这一项变成邻居 y 的加权平均ρ 的值被约束在可行区间里I - ρW也更容易保证可逆。两个最容易埋雷的细节一是对角线必须清零自己不能是自己的邻居二是距离矩阵如果不做行标准化谱半径可能超过 1后续求log|I - ρW|时直接得到 NaN。2.3 动态 SDM 与 SAR、SEM、静态 SDM 的选型关系模型的选取思路不是“越复杂越好”而是按照数据特征和理论问题一层层加上去。模型包含 y_t-1包含 Wy_t包含 WX典型估计策略回答什么问题SAR否是否ML邻居的 y 是否影响本地 ySEM否否否ML误差项是否存在空间相关静态 SDM否是是ML解释变量和被解释变量的双重溢出动态 SDM是是是QMLE / 空间 GMM短期效应与长期效应的区分我一般先跑静态 SDM用 LR 或 Wald 检验看能否退化为 SAR 或 SEM一旦发现被解释变量有明显的时间惯性比如上期污染物浓度对本期解释力很强再考虑动态 SDM。T 最好在 10 年以上T 太短时动态项和固定效应纠缠在一起估计结果极不稳定。3. 估计方法的取舍为什么 ML 在动态面板里会偏QMLE 与空间 GMM 怎么选3.1 内生性与 Nickell 偏误动态项的两个麻烦静态 SDM 常用的极大似然估计ML在动态面板里会遇到两个内生性问题。第一个来自 y_i,t-1个体固定效应 μ_i 进入所有期y_i,t-1 和 μ_i 必然相关这在经典动态面板里叫 Nickell 偏误。T 越小偏误越严重T 只有 10 到 20 期时τ 的估计会被明显向下拉。第二个问题来自 Wy_t空间滞后项是邻居当期的 y而邻居当期的 y 又受本地当期扰动影响Wy_t 和 ε_it 之间存在联立性也就是常说的“反射问题”。这两个问题叠在一起直接搬静态 ML 估计动态 SDM 的结果基本不能信。T 大到几十期时偏误才会逐渐消失但很多面板数据只有 10 到 15 年必须换估计策略。3.2 QMLE 偏误校正常见做法与算法骨架圈内最常见的做法是用带偏误校正的准极大似然估计QMLE。思路分三步先对面板做组内变换消掉个体固定效应再把空间滞后项和时间滞后项当作回归量用集中似然搜索 τ 和 ρ最后对组内变换导致的动态偏误做解析或自助法校正。很多现成的 Matlab 工具箱就是这么实现的。动手写自己的估计器之前至少先确认权重矩阵合格否则后面所有数值结果都是废的。我习惯在估计前跑一段矩阵检查% 检查行标准化后的空间权重矩阵 % W: N x N 空间权重矩阵要求对角线为零、行和为 1 maxEig max(abs(eig(full(W)))); if maxEig 0.9999 error(W 未正确行标准化谱半径接近或超过 1请检查对角线和对角线外的连接); end这段代码的逻辑是行标准化矩阵的最大特征值模理论上等于 1由于浮点计算允许一点点误差所以阈值设 0.9999。如果谱半径明显大于 1后面计算的log|I - ρW|在 ρ 较大时会直接变成无穷或 NaN。很多数据生成的流程跑着跑着出 NaN问题根本不在模型而在 W。3.3 什么时候换成空间系统 GMM如果 T 小于 15或者空间权重矩阵用的是经济距离这类与当期变量相关的矩阵QMLE 的偏误校正也不一定救得回来。此时可以考虑空间系统 GMM先做一阶差分消掉个体固定效应再用滞后两期及以上的水平项作为工具变量同时把 W*y 的空间滞后加入工具集。代价是工具变量数量容易爆炸。常见的翻车现场是 Sargan/Hansen 检验被拒绝原因往往不是模型错而是工具变量堆得太多或工具集里混入了当期内生项。空间系统 GMM 我一般只在 QMLE 结果明显异常、或者审稿人明确要求短面板稳健性时使用不会作为第一选择。4. 跑通动态面板空间杜宾模型的最小复现数据生成、权重矩阵与估计脚本4.1 准备一个可以自检的模拟数据DGP 先行不管读者是从别人手里拿到一个 .rar 压缩包还是攒了一套自己的面板数据第一步都是先让脚本能在一个“已知真实参数”的数据上跑通。Windows 下用 7-Zip 就能打开 rar 压缩包解压后先看里面有没有data、W、main三个部分如果没有就自己生成一份。我用一段可以完整运行的 Matlab/Octave 脚本生成动态 SDM 模拟数据。脚本里固定随机种子保证任何人复现的结果都一样rng(2024); % 固定随机种子保证结果可复现 N 60; T 10; % 60 个区域10 期面板 Wraw rand(N, N) 0.15; % 随机邻接关系平均每个区域约 9 个邻居 Wraw(1:N1:end) 0; % 对角线清零自己不是自己的邻居 W diag(1 ./ sum(Wraw, 2)) * Wraw; % 行标准化 Xmat randn(N, T); % 单个解释变量 mu randn(N, 1); % 个体固定效应 gam linspace(-0.5, 0.5, T); % 时期效应 tau 0.5; rho 0.3; % 真实参数 beta 1.0; theta 0.6; % 解释变量及其空间滞后系数 sig 0.3; % 扰动标准差 I speye(N); A I - rho * W; % 这一项必须可逆rho 要小于 1 Ymat zeros(N, T); for t 1:T xt Xmat(:, t); wxt W * xt; if t 1 rhs beta * xt theta * wxt mu gam(t) sig * randn(N, 1); else rhs tau * Ymat(:, t-1) beta * xt theta * wxt ... mu gam(t) sig * randn(N, 1); end Ymat(:, t) A \ rhs; end这段生成逻辑就是把动态 SDM 方程两边同时乘以(I - ρW)的逆把当期被解释变量显式解出来。t1 时没有上一期观测直接从静态分布启动之后每一期都同时受到上一期自己的影响、邻居当期的影响、解释变量和邻居解释变量的影响。参数说明tau0.5 表示时间惯性中等偏强rho0.3 表示空间溢出中等theta0.6 表示邻居的 x 对本地 y 有明显正向溢出。生成后建议先plot(Ymat(1:20,:))看一眼数据形态确认没有突变值。4.2 权重矩阵的行标准化与谱半径检查模拟数据里的 W 是随机生成的真实数据的 W 通常需要从地理坐标或行政区划构建。无论怎么构建估计前都按第 3 章的检查脚本跑一遍对角线是否为零、行和是否为 1、谱半径是否接近 1。如果发现某一行全为零说明存在孤立区域行标准化时会除零要提前处理成“以自身为邻居”或者删除该个体。% 检查是否有全零行 rowSum sum(Wraw, 2); if any(rowSum 0) warning(发现全零行请检查地理数据中是否有孤立区域); end全零行在邻接矩阵里很常见比如某个海岛县和任何陆地都不接壤。处理方式通常是两个要么用距离矩阵替代邻接矩阵要么把该区域从样本中剔除。强行保留会让行标准化那一步出 NaN后面的估计直接中断。4.3 QMLE 网格搜索脚本参数与初始值设置这里给出一版“集中似然 网格搜索”的完整脚本。它不是最前沿的偏误校正实现但把估计流程摆在了明面上对每个 τ 和 ρ先用外生变量回归得到 β、θ再计算残差和似然值最后取整个网格里似然最大的参数组合。真实研究中可以直接用现成的 QMLE 函数替换中间部分网格搜索的价值在于让人看清参数之间的耦合关系也方便排查为何优化器不收敛。% 输入Ymat(N x T)、Xmat(N x T)、W(N x N) % 输出tauHat、rhoHat、betaHat、thetaHat NT N * T; Yv Ymat(:); % 构造滞后项与空间滞后项 YlagMat [zeros(N, 1), Ymat(:, 1:end-1)]; WyMat W * Ymat; WXMat W * Xmat; % 组内去均值函数减去每个个体在 T 期上的均值 dm (z) z - kron(mean(reshape(z, N, T), 2), ones(T, 1)); Ydm dm(Yv); YlagDm dm(YlagMat(:)); WyDm dm(WyMat(:)); Xdm dm(Xmat(:)); WXDm dm(WXMat(:)); Z [Xdm, WXDm]; % 外生部分x 和 Wx % 参数网格先宽后密 tauGrid -0.5:0.05:0.9; rhoGrid -0.4:0.05:0.7; LL zeros(numel(tauGrid), numel(rhoGrid)); for a 1:numel(tauGrid) for b 1:numel(rhoGrid) tauVal tauGrid(a); rhoVal rhoGrid(b); % 构造剔除动态项后的因变量 ytmp Ydm - tauVal * YlagDm - rhoVal * WyDm; % 对剩余外生变量做 OLS得到 beta、theta bOLS Z \ ytmp; e ytmp - Z * bOLS; s2 (e * e) / (NT - size(Z, 2)); % 对数似然动态项贡献 空间雅可比项 logdet T * log(abs(det(I - rhoVal * W))); LL(a, b) -NT / 2 * (log(2 * pi * s2) 1) logdet; end end % 取网格最大值 [~, idxMax] max(LL(:)); [aIdx, bIdx] ind2sub(size(LL), idxMax); tauHat tauGrid(aIdx); rhoHat rhoGrid(bIdx); % 用最优 tau/rho 重新回归外生部分 ytmp Ydm - tauHat * YlagDm - rhoHat * WyDm; bOLS Z \ ytmp; betaHat bOLS(1); thetaHat bOLS(2); fprintf(tau_hat%.3f, rho_hat%.3f\nbeta_hat%.3f, theta_hat%.3f\n, ... tauHat, rhoHat, betaHat, thetaHat);代码里最关键的一行是ytmp Ydm - tauVal * YlagDm - rhoVal * WyDm它把时间滞后和空间滞后从因变量里“剔除”后剩余部分全部归因于外生变量。β 和 θ 的估计因此依赖于当前网格点的 τ 和 ρ这是集中似然的本质。网格范围要覆盖真实值的合理邻域动态面板 τ 一般落在 (-0.5, 0.9) 之间ρ 建议先参考静态 SDM 估计值再定。需要注意脚本保留了每个个体的第一期观测而严格研究中通常去掉第一期或者做初值处理否则滞后项的缺失值会影响 τ。4.4 从估计结果到短期/长期效应一张表说清楚空间杜宾模型里β 和 θ 不能直接当边际效应汇报。解释变量对 y 的影响会通过(I - ρW)^-1这个矩阵在网络里传播开所以要把直接效应、间接效应和总效应拆出来。单解释变量情况下短期效应矩阵是(I - ρW)^-1 * (βI θW)长期效应矩阵是((1-τ)I - ρW)^-1 * (βI θW)。% 从估计结果计算短期与长期效应 betaHat bOLS(1); thetaHat bOLS(2); Sshort (I - rhoHat * W) \ (betaHat * eye(N) thetaHat * W); Slong ((1 - tauHat) * I - rhoHat * W) \ (betaHat * eye(N) thetaHat * W); dirShort mean(diag(Sshort)); indShort mean(sum(Sshort, 2) - diag(Sshort)); dirLong mean(diag(Slong)); indLong mean(sum(Slong, 2) - diag(Slong)); fprintf(短期直接 %.3f间接 %.3f总效应 %.3f\n, ... dirShort, indShort, dirShort indShort); fprintf(长期直接 %.3f间接 %.3f总效应 %.3f\n, ... dirLong, indLong, dirLong indLong);这段代码的逻辑是短期效应矩阵只考虑当期空间反馈长期效应矩阵额外把时间累积算进去所以 τ 越大长期和短期的差距越大。汇报时可以整理成下面这种表审稿人最买账效应类型短期长期直接效应1.0231.214间接效应0.3870.562总效应1.4101.776注意不要把 βHat 填进“直接效应”一栏那是新手最容易犯的错。5. 动态 SDM 的五个常见坑现象、原因、解决一条条过5.1 似然值变成 NaN或者 ρ 被卡在网格边界现象网格搜索跑到一半log|I - ρW|输出 NaN或者最优 ρ 落在网格的最边缘比如 0.7 还想要更高。原因权重矩阵没有行标准化谱半径大于 1或者对角线没有清零自己和自己的空间滞后权重过强。解决先跑第 3.2 节的谱半径检查确认max(abs(eig(W)))在 1 附近。如果确实有个别行和不是 1重建矩阵不要手工缩放个别行直接统一行标准化。真实数据里如果存在全零行先处理孤立区域。5.2 滞后项串到别人的个体上τ 估计为负现象τ 的估计值和真实预期完全相反比如一个明显正相关的惯性过程估计出 -0.3并且 Wy 和 y_lag 的相关系数非常奇怪。原因面板数据从宽格式转长格式时个体排序没有按“个体-时间”排列滞后项[zeros(N,1), Y(:,1:T-1)]里的列顺序和当前行顺序错位前一期的数据跑到了别的个体上。解决进入估计前用sortrows或等价命令按(个体ID, 时间)排序排序后再构造YlagMat。我还会画前 20 个观测的连线图观察是否有断点跨过个体边界。5.3 把 β 系数直接当边际效应写进论文现象模型估计完正文说“解释变量每增加 1 单位被解释变量平均增加 β 个单位”。审稿人一眼看出问题。原因SDM 类模型通过(I - ρW)^-1把本地影响反馈到邻居再反馈回本地β 只是结构参数不是乘数结果。解决所有解释变量的汇报统一使用直接效应、间接效应和总效应不接受“先看系数再解释”。哪怕总效应方向和 β 相反也要以效应分解结果为准因为存在空间负反馈时 β 和总效应确实可能异号。5.4 动态项被当成“静态 SDM 加一个滞后变量”现象模型设定写了动态 SDM估计方法却填的是普通 ML且正文只讨论了 ρ 和 τ 的显著性没有讨论内生性处理。原因很多软件不能直接估计动态 SDM研究者就在静态模型里手动加y_t-1列用普通 ML 跑误以为加了滞后就是动态。解决务必在方法部分写明处理内生性的手段。如果是 QMLE要说明组内变换和偏误校正如果是空间系统 GMM要说明工具变量集。审稿人对“动态模型用静态估计”极其敏感这一条避不开。5.5 长期效应比短期效应大几个数量级现象长期间接效应是短期的几十倍甚至几百倍表格完全没法解释。原因(1 - τ)和ρ的关系出问题。长期效应矩阵((1-τ)I - ρW)接近奇异比如 τ0.85、ρ0.3 时最小特征值接近零逆矩阵被放大。解决估计后打印min(eig(full((1-tauHat)*I - rhoHat*W)))如果小于 1e-6说明模型接近非平稳。这时先看数据是否需要差分或者考虑把解释变量换成滞后项重新估计。不要强行解释一个条件数爆炸的矩阵。6. 汇报动态 SDM 结果的最后一步长期效应别只“乘个系数”要做矩阵稳健性6.1 从矩阵稳定性理解长期效应长期效应不等于短期效应乘以一个折扣。短期效应矩阵是(I - ρW)^-1 * (βI θW)它把当期空间传播算清楚了长期效应把时间累积也叠上去变成((1-τ)I - ρW)^-1 * (βI θW)。也就是说τ 通过改变I - ρW附近矩阵的条件数来影响长期结果而不是简单地把短期结果乘上1/(1-τ)。因此回归后的第一件事不是看系数表而是检查长期效应矩阵是否稳定。我习惯把下面这段判断放进所有动态 SDM 脚本里eigLong eig(full((1 - tauHat) * I - rhoHat * W)); if min(real(eigLong)) 1e-8 warning(长期效应矩阵接近奇异模型可能不平稳请检查 tau 和 rho 的组合); end这行命令比任何显著性检验都更能提前暴露问题。τ 和 ρ 单独看都显著合在一起却可能构成非平稳过程这类结果即使跑出来了也不能写进论文。6.2 权重矩阵敏感性换三种 W 看效应方向是否一致空间计量最容易被质疑的就是“换个权重矩阵结论就没了吗”。论文里一般要求做敏感性检验常见做法是准备三套矩阵邻接矩阵、距离倒数矩阵、K 近邻矩阵分别跑同一套估计脚本并输出三行效应。判断标准不是数值完全一致而是直接效应方向和总效应方向不要反转。如果邻接矩阵下溢出为正、距离矩阵下溢出为负说明结论对 W 过于敏感多半是模型设定漏了重要变量。我早期跑动态 SDM 时最喜欢盯着 τ 和 ρ 的显著性结果被审稿人一句“长期效应矩阵条件数太大”打回。后来养成一个习惯每次估计完先打印三个数——W 的谱半径、长期效应矩阵最小特征值、长期短期效应比值三个数都正常再开始写解释。希望这个习惯也能帮到你。本文还有配套的精品资源点击获取