简介动态面板空间杜宾模型是处理时空依赖与个体异质性的重要方法围绕该方法整理的实现资料包面向空间计量经济学研究者与高年级统计、经济类学生用于解决传统静态空间模型无法刻画时间动态效应的问题。压缩包内共19个文件以12个m脚本为主覆盖模型估计、空间权重矩阵构建、结果输出等功能模块另含4个Excel数据表格、1个PDF说明文档及少量临时文件PDF文档详细讲解极大似然估计原理数据表格可支撑实证案例复现整体仅337KB。资源结合高技术产业生产率与集聚等实际案例演示从数据准备、模型设定到运行与结果解读的完整路径有助于理解动态空间杜宾模型与静态模型的区别并快速上手二次开发。已有1948人学习参考适合具备基础空间计量知识、希望扩展动态面板建模能力的科研人员与高年级学生。1. 动态面板空间杜宾模型论文里那句“考虑动态调整”究竟让你加什么做空间面板回归的人大概都经历过这个场面静态空间杜宾模型跑完系数显著、效应分解也漂亮审稿人一句“该模型未考虑被解释变量的时间滞后与动态空间溢出”直接把结论打回。你重新翻文献发现近五年区域经济、能源环境、房地产领域的顶刊论文几乎都把这个“静态”换成了动态面板空间杜宾模型。标题里的“动态面板”指的就是在被解释变量的时间滞后项之外再引入空间滞后被解释变量的滞后项也叫时空滞后项让空间溢出效应有一个调整过程而不是瞬时完成。这套模型解决的问题非常具体上一期的环境污染、房价、经济增长会不会影响本期相邻地区的上一期表现又会不会通过空间传导影响本地区静态 SDM 假设所有调整瞬间完成动态版本放松了这个假设同时在效应分解上给出“短期”和“长期”两套直接效应与间接效应这正好是实证论文最需要的输出。适合的人群也明确手头有省份、城市或企业面板数据想用空间计量发论文或者被导师要求把静态结果升级成动态结果的人。下面这套流程按 Stata 和 Matlab 两条路展开都是实际能跑通的做法。2. 先把模型写对再动手动态空间杜宾的方程、滞后项与三种效应口径2.1 静态 SDM 到动态 SDM多出来的两个滞后项到底代表什么静态空间杜宾模型的标准写法是y_it ρΣ_j w_ij y_jt x_it β Σ_j w_ij x_jt θ μ_i ξ_t ε_it这里的 w_ij 是空间权重矩阵 W 第 i 行第 j 列元素ρ 是空间自回归系数x_it β 是本地区解释变量的直接影响Σ_j w_ij x_jt θ 是邻居解释变量带来的空间溢出。静态版本的经济含义是本地区 y 同时受当期邻居 y 和当期邻居 x 的影响且这种影响瞬间到达。动态面板空间杜宾模型在等式右边加上两项φ y_i,t-1 η Σ_j w_ij y_j,t-1第一项是时间滞后项反映状态依赖上一期的本地区 y 会对本期产生惯性作用比如污染排放的累积效应、房价的黏性。第二项是时空滞后项反映邻居上一期 y 经由空间关系传导到本地区本期 y 的动态溢出。注意时空滞后和静态模型里的空间滞后 Σ_j w_ij y_jt 有本质区别前者滞后一期后者是同期。这两个项引入后模型估计结果里会多出两个系数φ 是时间依赖强度η 是时空依赖强度。判断模型合不合理先看这两项是否显著。如果都不显著说明数据里没有明显动态效应强行做动态反而损失效率。我一般先跑一版静态 SDM 做基准再跑动态版用似然比检验或直接看 AIC 变化决定最终报告哪个设定。2.2 空间权重矩阵 W 的设定与行标准化所有结果的地基动态空间杜宾模型的全部空间行为都通过 W 来表达W 设错了ρ、η、θ 全是错的。空间权重矩阵分三类最常用邻接矩阵地理相邻取 1否则取 0、距离倒数矩阵w_ij 1/d_ijd_ij 为地区中心距离、经济距离矩阵用地区人均 GDP 差值的倒数构造。邻接矩阵最稳妥距离倒数矩阵适合解释力随距离衰减的场景经济距离矩阵审稿人经常质疑内生性——因为经济量本身受 y 影响。拿到或构造好 W 之后第一步永远是行标准化也就是让 W 每一行加起来等于 1。行标准化后的 W空间权重 w_ij 的含义从“是否相邻”变成“邻居影响的相对占比”同时保证 (I_N − ρW) 在合理 ρ 范围内可逆。我在 Stata 里做行标准化的常见做法是先用 Excel 或 Python 把原始权重表整理成 n×n 的 CSV然后在 Stata 里用 import delimited 读入再用 mkmat 转成矩阵最后逐行除以行和。提示权重矩阵行标准化必须在任何估计之前完成。很多人导入 0-1 邻接矩阵后直接喂给命令结果也能跑但效应分解数量级完全不对后面排查时很难发现是这里的问题。2.3 直接效应、间接效应与长短期分解审稿人最常追问的表格怎么来空间杜宾模型的回归系数不能直接解读为边际效应因为 W x 项的系数 θ 和空间溢出ρ会联动影响 y。正确的解读方式是做偏微分分解LeSage 和 Pace 给出过标准做法。以第 k 个解释变量为例静态 SDM 的效应矩阵是M_k (I_N − ρW)^{-1} (β_k I_N θ_k W)这个矩阵的对角线均值是平均直接效应含义是本地区第 k 个变量 x 对本地区 y 的平均影响非对角线元素反映一个地区 x 变化对另一个地区 y 的影响所有非对角线元素求平均得到平均间接效应也就是空间溢出效应。直接效应加间接效应等于总效应。加入时间滞后项 φ 之后长期效应需要在静态效应矩阵基础上乘一个调整因子 1/(1−φ)。这个逻辑很直观y_it 的变化会通过 φ y_i,t-1 递推到下一期累积起来相当于放大了 1/(1−φ) 倍。如果同时加入了时空滞后项 η长期效应矩阵变成M_long,k ( (1−φ)I_N − (ρη)W )^{-1} (β_k I_N θ_k W)不能简单乘 1/(1−φ) 了事因为 η 和 ρ 合在一起出现在矩阵求逆里。这也是我见过最多的翻车点很多人不管 η 是否显著一律乘 1/(1−φ)导致长期溢出效应高估或低估。记住一个硬性条件模型必须满足 (1−φ) 与 ρη 的组合落在特征值约束区间内否则长期矩阵发散结果出现 NaN 或极端大数。3. 用 Stata 跑通动态面板空间杜宾xsmle 命令与一份可直接套用的回归流程3.1 数据准备平衡面板、权重矩阵 CSV、变量检查Stata 里做动态空间杜宾估计xsmle 是应用最广的外部命令支持空间滞后模型、空间误差模型和空间杜宾模型并且 dlag 选项专门处理动态设定。第一步安装ssc install xsmle, replace数据准备阶段有三个硬要求。第一面板必须是强平衡面板每个个体有相同的期数不允许缺年第二权重矩阵 W 必须是 n×n 的 Stata 矩阵对象不能是变量第三被解释变量不能有缺失值否则滞后项构造会传染。权重矩阵导入的完整流程这样写假设有 30 个地区、权重表存在 W_30.csv 里import delimited W_30.csv, clear mkmat v1-v30, matrix(W) save W_matrix, replacemkmat 后面的 v1-v30 来自 CSV 导入后自动命名的变量顺序必须和面板数据里的个体 id 顺序一致。我踩过一次坑权重表按省份名称排序面板数据按 id 排序两个顺序对不上估计结果里空间系数虽然显著但方向完全相反。所以导入后第一件事是核对 W 的第一行和面板里 id1 对应的是不是同一个地区。3.2 xsmle 估计dlag(1)、dlag(2) 与双向固定效应的完整命令数据整理成 xtset 格式后用 xsmle 跑回归。先看一个完整命令use panel_data.dta, clear xtset id year global xlist x1 x2 x3 * 动态SDM仅含时间滞后项 xsmle y $xlist, wmat(W) model(sdm) dlag(1) fe type(both) nolog * 动态SDM含时间滞后和时空滞后项 xsmle y $xlist, wmat(W) model(sdm) dlag(2) fe type(both) nolog第一段命令里的 dlag(1) 表示只加入 φ y_i,t-1dlag(2) 表示同时加入 φ y_i,t-1 和 ηΣ_j w_ij y_j,t-1对应 2.1 里的完整动态模型。fe 表示固定效应type(both) 指定个体和时间双向固定效应这是空间面板论文里的默认设定因为不控制时间固定效应时宏观冲击会被错误归因到空间相关性上。wmat(W) 接收 Stata 内存里的矩阵 W。模型选择 model(sdm) 是因为 SDM 是嵌套了空间滞后模型和空间误差模型的更一般形式如果检验发现空间误差项占主导再退回 model(sem)。nolog 用来隐藏迭代过程输出更清爽。回归完成后保存估计结果方便后续比较est store dyn_sdm_dlag23.3 读懂 xsmle 输出从 rho、滞后项系数到短期/长期效应表xsmle 运行后结果表上半部分是解释变量系数下半部分有几行特殊内容需要单独看。以 dlag(2) 的估计为例rho本期空间自回归系数衡量同期邻居 y 的影响y_lag时间滞后项系数 φwy_lag时空滞后项系数 ηWx 开头的行解释变量空间滞后的系数 θ如果 y_lag 和 wy_lag 都显著为正说明被解释变量存在明显惯性和空间传导惯性。此时需要进一步看短期与长期效应分解。xsmle 在动态设定下会自动报告 short-run 和 long-run 两套直接效应与间接效应表格底部会额外输出这些结果不用额外命令。从结果整理成论文表格时我通常手动把短期直接效应、短期间接效应、长期直接效应、长期间接效应拼成四列每列下面放 t 值或置信区间。长期效应和短期效应的差异大小本身就是一个可写进论文的点如果长期间接效应远大于短期间接效应说明空间溢出在时间维度上缓慢累积政策评估如果只报短期值会严重低估影响。3.4 显著性与稳健性固定效应和随机效应的选择审稿人经常会质疑固定效应模型里个体效应和时间效应的设定。常规做法是跑一版 fe、跑一版 re然后用 Hausman 检验给结论。xsmle 支持随机效应估计xsmle y $xlist, wmat(W) model(sdm) dlag(2) re nolog est store dyn_sdm_re hausman dyn_sdm_dlag2 dyn_sdm_re这里要注意 hausman 命令要求两个模型估计同样的解释变量集合并且权重矩阵一致。Hausman 检验的原假设是随机效应与解释变量不相关如果 p 值小于 0.05拒绝原假设报告固定效应否则报告随机效应。实际经验里空间面板的随机效应估计在 xsmle 里偶尔会遇到收敛警告如果 re 版本不收敛通常是在解释变量里加入了地理属性变量导致与个体效应多重共线删除该变量即可。4. 换 Matlab 做动态空间杜宾jplv7 流程与效应分解的一个可复现模板4.1 为什么还要留一手 Matlab动态模型里 xsmle 的局限xsmle 能覆盖大多数情况但有一个短板它对长期效应的标准误计算依赖 delta 方法结果表不直接给出长期效应的诊断信息且 dlag(2) 的稳定性约束处理得比较隐讳。另一个现实问题是很多审稿人和看论文的复现者更熟悉 Elhorst 和 LeSage 那套 Matlab 代码尤其 jplv7 工具箱在国内学术圈流传很广。手里备一套 Matlab 流程一方面能交叉验证 Stata 结果另一方面当 xsmle 因面板不平衡或收敛问题报错时有一个退路。jplv7 是 LeSage 公开的空间计量工具箱里面包含大量面板估计函数和辅助工具函数最常用的是 sar_panel 系列。动态空间杜宾模型不直接对应某个单一函数常见做法是把滞后项构造好再交给静态面板空间估计器处理。这套流程的核心不是某一段黑匣子代码而是下面三个模块化的步骤。4.2 数据排列与滞后项构造N×T 堆叠格式下的两个自定义函数jplv7 的风格是数据按“个体堆叠”排列前 T 行是第一个个体的全部时期接着 T 行是第二个个体。假设有 N 个地区、T12 期y 的长度是 NT×1。构造时间滞后项不能直接用 MATLAB 内置 lag 命令它会跨个体滚动必须按个体内滞后。自定义函数如下function ylag lag_byid(y, T) % y 是 NT×1T 是每个个体的期数 % 输出 ylag(i,t) y(i,t-1)每个个体第一期补 0 N length(y) / T; ylag zeros(size(y)); for ii 1:N idx (ii-1)*T 1 : ii*T; ylag(idx(2:end)) y(idx(1:end-1)); end end这个函数保证每个个体内部错位一期不跨个体。时间滞后项构造好后空间滞后和时空滞后通过 Kronecker 积实现Wbig kron(W, eye(T)); % 空间权重扩展到 NT×NT wy Wbig * y; % 同期空间滞后 wylag Wbig * ylag; % 时空滞后 W*y(t-1)kron(W, eye(T)) 的含义是把 W 里的每个元素扩成 T×T 的单位矩阵块使得乘上去之后每个地区本期 y 对应同一时期邻居的 y。参数说明W 必须先用行标准化eye(T) 保证时间维度上不混叠。这一步是 Stata 里 xsmle 自动完成的在 Matlab 里必须手动处理。4.3 估计主程序模型矩阵拼装、normw 行标准化与 ML 估计入口估计前的准备工作是把动态项并入解释变量矩阵。以两个核心解释变量为例% 数据准备 N 30; T 12; K 2; y data(:, 1); % NT×1 x data(:, 2:3); % NT×K % 权重矩阵行标准化 W normw(W); % jplv7 自带行和归一为 1 % 或手动标准化W W ./ sum(W, 2); % 构建设定项 ylag lag_byid(y, T); Wbig kron(W, eye(T)); wy Wbig * y; wylag Wbig * ylag; Wx Wbig * x; % 拼装为预定义变量矩阵 Xall [x Wx ylag wylag]; % 调用 jplv7 的空间固定效应面板估计器 info.lflag 0; % 完整似然而不是近似 info.model 1; % 固定效应 results sar_panel_FE(y, Xall, W, T, info);代码里 normw 是 jplv7 自带的行标准化函数也可以手动逐行除以行和。info.lflag0 指定计算精确对数似然样本太大时建议改为 lflag1 启用近似速度差好几倍。sar_panel_FE 内部会对 Xall 做面板变换把固定效应消去后再做 ML 估计。有一个必须直说的注意点直接把 ylag 和 wylag 当普通解释变量放进 ML 估计器等于把动态项视为预先确定变量严格来说这是简化处理。严谨的动态空间面板需要 GMM 或基于偏差校正的 ML 方法。但对复现和初筛来说这个简化结果足以判断动态效应方向、显著性以及长期效应量级正式投稿前再移植到完整的动态 ML 程序里。4.4 长短效应分解把公式翻译成矩阵运算估计完成后从 results 里提取 ρ、β、θ、φ、η再做效应分解。核心代码如下rho results.rho; beta results.beta(1:K); % 前 K 个为 x 的系数 theta results.beta(K1:2*K); % 中间 K 个为 Wx 的系数 phi results.beta(2*K1); % 时间滞后系数 eta results.beta(2*K2); % 时空滞后系数 S inv(eye(N) - rho * W); I_N eye(N); direct_sr zeros(K, 1); indirect_sr zeros(K, 1); direct_lr zeros(K, 1); indirect_lr zeros(K, 1); for k 1:K Mk_sr S * (beta(k) * I_N theta(k) * W); direct_sr(k) trace(Mk_sr) / N; indirect_sr(k) (sum(Mk_sr(:)) - trace(Mk_sr)) / N; % 长期效应需解 (I - phi*I - (rhoeta)*W) 的逆 LongInv inv((1 - phi) * I_N - (rho eta) * W); Mk_lr LongInv * (beta(k) * I_N theta(k) * W); direct_lr(k) trace(Mk_lr) / N; indirect_lr(k) (sum(Mk_lr(:)) - trace(Mk_lr)) / N; end短期效应矩阵 S 的推导来自 (I_N − ρW)^{-1}。长期效应的 LongInv 对应 2.3 里的公式注意必须代入 φ 和 η 的估计值不能用静态矩阵。sum(Mk(:)) 是矩阵所有元素之和取平均后得到总效应减去迹平均就是间接效应。 run 完这段用表格形式把短期/长期直接和间接效应对齐到每个变量和 Stata 的 xsmle 输出对比如果两边的间接效应方向或量级差异超过 20%优先怀疑权重矩阵标准化或数据排列顺序问题。5. 动态面板空间杜宾的五个高频坑从 rar 解压到结果对不上的排查手册5.1 rar 压缩包解压报错、提示密码先看注释和 Readme别急着找移除工具很多读者拿到的是类似“动态面板空间杜宾模型.rar”这种分享包包内文件名还带着 caughtuk3、v2_final 这类个人标记。现象是用 WinRAR 或 7-Zip 解压时提示文件头损坏、或者弹密码输入框再或者在中文系统下解压出乱码文件名。原因通常是三选一下载不完整导致压缩包末尾截断作者打包时设了密码但说明写在压缩包注释里文件名用了非 UTF-8 编码系统解码错乱。解决步骤先别急着找什么 rar密码移除、recovery toolbox 破解版之类的工具那些对加密包的恢复基本无效而且这类工具是捆绑病毒的重灾区。正确做法是用 7-Zip 打开压缩包先看右侧注释栏有没有作者留下的密码或说明再试一次完整重新下载并比对文件大小是否和分享页一致文件名乱码时在 7-Zip 里通过“工具-选项-编码”切换查看编码。密码移除软件不仅帮不上忙还可能让论文数据和电脑一起翻车。如果包内只有代码没有数据问题不大如果数据也在包里且解压失败老实联系分享者要密码最靠谱。5.2 权重矩阵忘记做行标准化结果表面上能跑效应分解一算就翻车现象xsmle 和 Matlab 端都能正常收敛ρ 也显著但间接效应大得离谱或者直接效应和总效应符号相反。原因把原始 0-1 邻接矩阵直接传入行和不等于 1导致 (I_N − ρW)^{-1} 里空间乘子的缩放尺度错误。xsmle 的帮助文档明确要求权重矩阵行标准化但它不会主动检查你传什么它用什么。手动改法在 Stata 里构建矩阵后跑一个简单的 Mata 循环或直接在 Excel/CSV 阶段完成标准化。在 CSV 阶段最省事的方法是每行数值除以该行总和形成一个新表再导入。提示标准化的本质是让权重矩阵的特征值上限变得可控。验证方法很简单对标准化后的矩阵计算最大特征值理论上接近 1ρ 的有效取值范围是 1/λ_min 到 1/λ_max 之间超出就会出现奇异矩阵警告。5.3 xsmle 报 strongly balanced 错误非平衡面板怎么处理现象数据里有几个地区起始年份晚一年或者中间缺了一期xtdes 显示非平衡xsmle 直接报错。原因动态模型要构造 y_i,t-1缺期会产生无法估值的缺口xsmle 要求必填的平衡面板数据里有一行的 year 与 id 组合缺失都会触发这个错误。解决先用 xtset 和 xtdes 检查缺漏再用 xtbalance 截尾xtset id year xtbalance, range(2005 2020) save panel_balanced.dta, replacextbalance 会把样本限定在共同的时间区间内有缺失年的个体整体剔除。注意fillin 补齐缺失年份再插值在空间动态模型里是下策因为插值会人为制造空间相关审稿人一旦追问数据构造过程很难解释。宁可损失一点样本量也不要制造假观测。5.4 Matlab 跑出 NaN、复数或发散滞后期太长、初始值丢失和 (I−ρW) 奇异性排查现象sar_panel_FE 运行后 beta 出现 NaN 或复数值有时提示矩阵接近奇异或尺度太大。原因有几种按出现频率排序数据里 y 或 x 含 NaN滞后函数 lag_byid 会把 NaN 逐期向后传递导致一半样本被污染权重矩阵没有标准化ρ 搜索时 (I − ρW) 接近奇异导致对数似然函数在迭代边界上失效T 和 N 的比例失调比如 N10、T2动态项太多导致识别不足。解决顺序先检查数据里有没有 NaN再检查 eig(W) 的最大特征值最后限制 ρ 的搜索区间。jplv7 的 sar_panel_FE 支持 info 里的 rmin 和 rmax 参数手动设置info.rmin -1.2 / max(abs(eig(W))); info.rmax 1.2 / max(abs(eig(W)));这个设置给 ρ 留一个不超过谱边界的区间避免迭代过程中越过奇异点。复数结果基本可以断定是迭代步长跨过了不可逆区域调小 rmax 是最快解法。5.5 Stata 与 Matlab 结果“差不多但不一样”权重矩阵读入方向与数据排列的差异现象同一个数据和权重矩阵xsmle 和 sar_panel_FE 的结果方向一致但数值对不上ρ 差 0.1、间接效应差 30%。原因极大概率是权重矩阵排列方向W(i,j) 在数学定义里表示 j 对 i 的影响但导入时如果行和列被转置空间滞后的方向就反了。另一个常见原因是数据堆叠顺序不一致Stata 的 xtset 是 id 优先排序而 Matlab 脚本里如果读入数据时按时间优先排列kron(W, eye(T)) 就会错误地把权重乘到错误的时间段上。解决写一段核对代码在两种环境里同时打印 W 第一行前五个元素和 y 第一个个体的前三个观测值确认矩阵和数据结构完全一致。再各跑一个简化模型只保留一个解释变量比较 ρ 和 β 是否一致逐步加项排查差异来源。这个方法朴素但效率最高比反复调参数靠谱得多。6. 权重矩阵跑三套、长短期效应一起报动态空间杜宾结果的最后一道自检动态空间杜宾模型在审稿阶段最容易受到的攻击是“结果对权重矩阵选择敏感”。我的应对习惯是主回归用 Queen 邻接矩阵稳健性检验用距离倒数矩阵和经济距离矩阵各跑一遍并且保证动态项设定不变。三套矩阵下ρ 的符号、动态项 φ 和 η 的显著性、长期间接效应的方向必须保持一致只要有一套出现符号翻转哪怕主回归再漂亮也要回去检查是权重构造问题还是模型设定问题。具体操作上我会在脚本里把三套权重矩阵的构建和标准化写成同一个函数确保除了矩阵本身以外所有估计参数完全一致避免复制粘贴导致的不必要误差。表格里报告结果时把三套矩阵下的长期直接效应和长期间接效应并列排放旁边加一行特征值范围说明让审稿人一眼看到模型在参数空间内是稳定的。另一个自检技巧是把结果与静态 SDM 做对比如果动态项不显著长期效应和静态效应应该非常接近如果动态项显著长期直接效应通常大于静态直接效应。要是出现长期效应小于静态效应的反常情况回头检查 2.3 里的长期效应公式是否用对了。我自己保留的习惯是把数据、权重矩阵、Stata do 文件和 Matlab m 文件放在同一个目录权重矩阵文件名里带上版本号比如 W_adjacent_v2.csv重新跑复现时先核对文件版本而不是急着跑回归。这个习惯救过我很多次。希望帮到你。本文还有配套的精品资源点击获取