动态面板空间杜宾模型实战:从方程设定到Stata/MATLAB实现
2026/9/24 23:58:43 网站建设 项目流程

简介:动态面板空间杜宾模型是处理时空依赖与个体异质性的重要方法,围绕该方法整理的实现资料包面向空间计量经济学研究者与高年级统计、经济类学生,用于解决传统静态空间模型无法刻画时间动态效应的问题。压缩包内共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_ij,d_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, replace

mkmat 后面的 v1-v30 来自 CSV 导入后自动命名的变量,顺序必须和面板数据里的个体 id 顺序一致。我踩过一次坑:权重表按省份名称排序,面板数据按 id 排序,两个顺序对不上,估计结果里空间系数虽然显著但方向完全相反。所以导入后第一件事是核对 W 的第一行和面板里 id=1 对应的是不是同一个地区。

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-1;dlag(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_dlag2

3.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 个地区、T=12 期,y 的长度是 NT×1。构造时间滞后项不能直接用 MATLAB 内置 lag 命令,它会跨个体滚动,必须按个体内滞后。自定义函数如下:

function ylag = lag_byid(y, T) % y 是 NT×1,T 是每个个体的期数 % 输出 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.lflag=0 指定计算精确对数似然,样本太大时建议改为 lflag=1 启用近似,速度差好几倍。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(K+1:2*K); % 中间 K 个为 Wx 的系数 phi = results.beta(2*K+1); % 时间滞后系数 eta = results.beta(2*K+2); % 时空滞后系数 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 - (rho+eta)*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", replace

xtbalance 会把样本限定在共同的时间区间内,有缺失年的个体整体剔除。注意:fillin 补齐缺失年份再插值,在空间动态模型里是下策,因为插值会人为制造空间相关,审稿人一旦追问数据构造过程很难解释。宁可损失一点样本量,也不要制造假观测。

5.4 Matlab 跑出 NaN、复数或发散:滞后期太长、初始值丢失和 (I−ρW) 奇异性排查

现象:sar_panel_FE 运行后 beta 出现 NaN 或复数值,有时提示矩阵接近奇异或尺度太大。原因有几种,按出现频率排序:数据里 y 或 x 含 NaN,滞后函数 lag_byid 会把 NaN 逐期向后传递,导致一半样本被污染;权重矩阵没有标准化,ρ 搜索时 (I − ρW) 接近奇异导致对数似然函数在迭代边界上失效;T 和 N 的比例失调,比如 N=10、T=2,动态项太多导致识别不足。

解决顺序:先检查数据里有没有 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,重新跑复现时先核对文件版本而不是急着跑回归。这个习惯救过我很多次。希望帮到你。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询