简介：这份资源面向空间计量经济学、区域经济与社会科学领域的研究者，聚焦动态面板空间杜宾模型的实现与应用。相比传统静态空间模型，动态设定同时纳入时间滞后与空间滞后效应，能更稳健地刻画技术溢出、产业集聚等长期动态关系。资源内含12个Matlab脚本文件，覆盖模型估计、空间权重构建、结果输出及配套示例；4个Excel文件用于存放研究数据，并附有相关PDF说明文档，可帮助用户理解估计思路与编程细节。压缩包共19个文件，整体仅337KB，轻量便携，便于直接下载后运行示例并替换为自有数据开展实证分析。目前已有1949人学习下载，适合具备一定空间计量基础、希望从静态模型转向动态模型或快速上手Matlab代码的学者与学生使用。

1. 动态面板空间杜宾模型：静态模型看不出的时空交互，才是这份资源的重点

做高技术产业生产率与集聚的空间面板分析时，很多人第一反应是静态空间杜宾模型（SDM）。但静态模型只捕捉当期的空间溢出，而产业集聚、技术扩散天然带滞后性——企业今天的区位选择，往往受邻区上一期生产率的牵引，这种时空交互静态模型看不见。这份“动态面板空间杜宾模型”资源包正好补上这块：从 note-sdpd-mle.pdf 推导笔记，到 sar_jihai、SARAR 等可直接跑的 MATLAB 函数，再到高技术产业面板数据样本，一条龙齐备。论文党、区域经济研究者只要数据是“多个地区×多个年份”，就能从静态平滑切到动态。这篇文按我的复现习惯，把建模、权重处理、参数设置到避坑的完整过程拆开讲。

2. 先看懂模型再跑代码：动态SDM的设定与MLE原理

2.1 动态项来自哪里：从静态SDM到动态SDM的方程推演

静态空间杜宾模型的标准写法是：

y = ρWy + Xβ + WXθ + ε

其中 ρ 是同期空间溢出系数，WXθ 是解释变量的空间滞后项，用来捕捉“邻居的解释变量对本地被解释变量”的影响。它比 SAR 模型多一个 WX 项，比 SEM 模型更适合单独识别溢出效应，所以这几年实证论文里用得很多。

但静态 SDM 有一个隐含假设：所有交互都发生在同一时期。这对经济现实来说过于勉强。拿高技术产业来说，企业学习邻居的先进技术需要时间，政策从出台到影响企业决策也有时滞；同时上一期本地的高生产率本身也会在下一期继续显现效应。这两件事静态模型都装不进去。

动态 SDM 在等式右侧补了两项：

y_t = τy_{t-1} + ρWy_t + ηWy_{t-1} + X_tβ + WX_tθ + μ_i + ε_t

与静态版本逐项对比：

• τ：上一期本地区被解释变量对本期的影响，也就是时间惯性项。
• η：上一期邻居被解释变量对本期本地区的影响，这是时空扩散项，也是动态 SDM 区别于静态模型的核心所在。
• ρ：本期邻居对本地区的影响，同期空间溢出。
• β 和 θ 沿用静态 SDM 的含义，分别是解释变量的直接效应和空间溢出效应。

当 τ 和 η 同时为 0，模型退化为静态 SDM；当 η 为 0，退化为动态 SAR；当 ρ 和 η 都为 0，退化为普通面板模型的动态回归。所以动态 SDM 是嵌套了静态 SDM 的更一般设定，也正因为嵌套，可以直接用似然比或 LM 检验来判断是否需要上动态版本。

我一般怎么决定要不要跑动态？看两个信号。第一，被解释变量的一阶自相关显著大于 0.8，说明时间惯性很强；第二，静态 SDM 回归后残差仍有空间自相关，但加 WX 项也没消除。这两个信号只要中一个，静态结果很容易被质疑遗漏时空滞后变量导致内生性，这时就该上动态。

2.2 似然函数与空间雅可比项：为什么MLE不是简单OLS

空间模型的核心难点在于 ρWy_t 与误差项相关，直接用 OLS 得到的 ρ 和 β 都有偏。所以标准做法是极大似然估计（MLE）。面板形式下，对数似然函数可以写成：

lnL = -NT/2 × ln(2πσ²) - T × ln|I_N - ρW| - 1/(2σ²) × e′e

其中 e 表示经过动态项、空间项和解释变量拟合之后的残差。式子里的 ln|I_N - ρW| 就是空间雅可比项，它把空间自回归系数 ρ 的估计拉回可行域内。没有这一项，估计出来 ρ 很可能落在 [-1, 1] 之外，经济含义直接失效。

问题在于，当 N 是省份数量 31 的时候，直接算矩阵行列式还能忍受；一旦换成城市级样本，N 上百甚至几百，每次优化迭代都算一次行列式，速度就拖垮了。常见做法是先对 W 做一次特征值分解：

ln|I_N - ρW| = Σᵢ ln(1 - ρλᵢ)

其中 λᵢ 是 W 的特征值。这样特征值只算一次，后续每次迭代只要做 N 次对数求和即可，优化速度提升一个量级。这个资源里的 f2_sar_jihai.m 和 f2_sar_jihai_time.m，本质上就是把给定 ρ 时的似然值或残差平方和封装成函数，供 fminbnd 这类一维搜索器迭代。

这类目标函数通常采用集中似然（concentrated likelihood）的思路：内层把 β 和 σ² 用 OLS 显式解出，外层只对 ρ 做一维搜索。这样把多维优化降成一维，稳定性高出不少。f2 前缀基本可以理解为第二版目标函数，比第一版更收敛。

提示：MLE 的搜索区间必须满足 ρ ∈ (1/λ_min, 1/λ_max)。如果 W 没做行标准化，λ_max 可能远超 1，搜索区间设置不合适时结果会直接卡在边界上。这个在后面避坑章节会专门展开。

2.3 资源包文件清单：每个文件在流水线里的角色

把整个包的文件摊开看，可以分成五组：

从命名风格看，这套代码沿用了常用空间计量工具箱的函数组织方式，但把动态项和时变逻辑加了进去。跑之前先读 note-sdpd-mle.pdf 是个好习惯，这份笔记把动态 SDM 的 MLE 推导、空间雅可比处理都写清楚了，等于给后续调参打底。

3. 空间权重矩阵与数据准备：权重错了，估计系数全白搭

3.1 wrook.m与Rook邻接矩阵的生成逻辑

空间权重矩阵是空间模型的第一道关卡。所谓 Rook 邻接，就是两个区域存在共享边界（上、下、左、右）才记为 1，对角接触不算；Queen 邻接则连对角也算。Rook 更严格，用于省级面板时不容易造成过度连接。

wrook.m 的思路是把 N 个区域按网格排列，然后对每个格子找四方向邻居：

% wrook.m：把N个区域排成n行m列网格，生成Rook邻接 % 输入：n 网格行数，m 网格列数 % 输出：N*N 稀疏邻接矩阵，N = n*m function W = wrook(n, m) N = n * m; W = sparse(N, N); for i = 1:N [r, c] = ind2sub([n, m], i); % 一维编号转二维网格坐标 if r > 1, W(i, i - n) = 1; end % 上方邻接 if r < n, W(i, i + n) = 1; end % 下方邻接 if c > 1, W(i, i - 1) = 1; end % 左方邻接 if c < m, W(i, i + 1) = 1; end % 右方邻接 end end

逻辑说明：稀疏矩阵只存非零元素，N 上百时内存占用大幅降低。ind2sub 把一维编号转成二维坐标后，上下邻居相隔 n，左右邻居相隔 1，四个判断条件覆盖完整的 Rook 邻接。

但对实际论文数据来说，wrook.m 更多是教学或规则网格场景才适用。省级面板区域不是规则网格，直接用邻接表或经济距离矩阵更常见。对这个包而言，matrix.xlsx 和“高技术产业生产率与集聚W.xlsx”才是真正的权重来源，wrook.m 只能作为生成初始邻接的参考模板。

3.2 从matrix.xlsx读取并做行标准化：特征值体检

实际研究中 W 矩阵通常已经事先整理好，放在 Excel 里。读取后做行标准化是必须的一步：

% 从Excel读取空间权重矩阵 W_raw = xlsread('matrix.xlsx'); % 体检1：必须是方阵 assert(isequal(size(W_raw, 1), size(W_raw, 2)), '权重矩阵必须是方阵'); % 行标准化：每行元素除以行和 W = W_raw ./ sum(W_raw, 2); % 体检2：对角是否为0（存在自环会污染rho） if any(diag(W) ~= 0) warning('matrix.xlsx 对角存在非零值，请检查自环问题'); end % 体检3：行标准化后每行和为1 rowsum = sum(W, 2); assert(max(abs(rowsum - 1)) < 1e-10, '行标准化失败');

参数说明：行标准化后每行和为 1，最大特征值刚好是 1，此时 ρ 的可行域是 (1/λ_min, 1)。不标准化的话，λ_max 可能是几十，ρ 的估计范围直接被撑爆，数值优化很容易翻车。

读完 Excel 之后做特征值体检是我固定保留的步骤：

lambda = eig(W); fprintf('特征值范围：%.3f 到 %.3f\n', min(lambda), max(lambda));

如果 max(lambda) 明显偏离 1，回去检查 Excel 而不是在 MATLAB 里调参。最常见的低级错误是 Excel 打开状态下直接复制数据，导致 MATLLAB 读到被锁定的临时文件；这个资源包里的“~$副本高技术产业生产率与集聚W.xlsx”文件就是这样残留的，打包前没关 Excel。

3.3 面板数据堆叠与动态滞后变量的构造

data.xlsx 里是“高技术产业生产率与集聚”的面板数据。用之前先确认排列顺序：省份升序，同省内部按年份升序。这个顺序直接决定滞后变量 y_{t-1} 是否准确。

% 读取面板数据 raw = xlsread('data.xlsx'); prov = raw(:, 1); % 省份编码 year = raw(:, 2); % 年份 y = raw(:, 3); % 被解释变量 x = raw(:, 4:end); % 解释变量矩阵 % 排序检查与重排 panel = sortrows([prov, year, raw(:, 3:end)], [1, 2]); y = panel(:, 3); x = panel(:, 4:end); N = length(unique(prov)); % 地区数 T = length(unique(year)); % 期数

逻辑说明：sortrows 按第 1 列省份、第 2 列年份升序排列，确保同一地区的时间序列按年份连续排列。这样后面做 y_{t-1} 时，每一行的滞后值才对应上一年的真实观测，而不是随便某个样本。

接下来构造动态回归需要的滞后项。这一步经常被忽略，但直接决定模型能不能跑出来：

% 构造滞后变量：每个地区上一年的y ylag = zeros(N * T, 1); for i = 1:N idx = (i - 1) * T + (1:T); % 该地区所有年份的下标 ylag(idx) = [NaN; y(idx(1:T-1))]; % 滞后一期，第一年置NaN end % 空间滞后 wy = W * reshape(y, N, T); % 逐期计算Wy wylag = W * reshape(ylag, N, T);

这里 ylag 的第一年对应位置是 NaN，估计动态模型时必须把每个地区第一年的观测剔除，否则 NaN 会传导到似然函数里，结果变成整片 NaN。这个数据清洗细节，是跑通动态 SDM 的前提。

提示：如果面板是非平衡的，比如某些年份缺失，滞后变量的下标会错位。先在 Excel 里补齐空缺年份，或者删掉非平衡地区，再进入估计。

4. 核心估计函数调用与参数设置：从sar_jihai到sar_jihai_time

4.1 sar_jihai.m：静态空间回归的调用与输出

sar_jihai.m 估计的是静态 SAR，也就是 y = ρWy + Xβ + ε。它作为动态模型的热身和基准，承担“先看静态结果，再判断要不要加动态项”的作用。调用格式如下：

% 静态SAR估计 % 输入： % y 被解释变量（NT x 1） % W 空间权重矩阵（N x N），已行标准化 % X 解释变量（NT x K），不含常数项 % T 期数 % 输出： % results 结构体，包含rho、beta、tstat等 results = sar_jihai(y, W, x, T); % 输出回归结果 prt_reg(results); prt_sp(results);

参数说明：T 必须传入，因为面板模型的误差结构按 T 期堆叠，似然函数里要对 T 求幂。X 不需要包含常数项，常数会在估计内部自动处理，手动加常数项反而会导致完全共线性，输出 rank deficient。

prt_reg 打印的是通用回归表，prt_sp 打印空间模型特有的 rho、对数似然值、AIC、BIC 这些。两个函数配合，适合做静态 SDM 与动态 SDM 的模型比较。跑完静态后看一眼 rho 的显著性：如果 rho 不显著，可能根本不需要空间模型，直接面板回归就行。

4.2 sar_jihai_time.m与SARAR.m：时变系数和误差空间相关的稳健性

某些研究里 ρ 并不是常数。比如样本期内发生了重大区域政策调整，空间溢出强度很可能前后不同，这时用整个时期一个 ρ 会掩盖结构变化。sar_jihai_time.m 解决的就是这个问题：

% 时变空间自回归：每个时间截面估计一个 rho_t results = sar_jihai_time(y, W, x, T); % 输出各期空间自相关系数 fprintf('各期空间自相关系数：\n'); disp(results.rho_t);

逻辑说明：每个 t 期单独估计 ρ_t，相当于 T 个截面各自跑 SAR 再汇总。代价是 T 必须足够大，一般建议不低于 10；单期样本量太小，ρ_t 的方差会膨胀到失去解释意义。如果 T 只有 5 年左右，不建议用这个函数，直接用固定 ρ 的动态 SDM 更稳。

SARAR.m 则是另一种稳健性检验思路。SARAR 模型在 SAR 基础上再加一个空间误差项的回归：

y = ρWy + Xβ + u，u = λWu + ε

它适用于被解释变量存在空间溢出、同时不可观测因素也存在空间相关的场景。如果跑完动态 SDM 后残差仍显示空间相关，可以用 SARAR.m 做一个对比版本：

% SARAR估计 results_sarar = SARAR(y, W, x, T); prt_sp(results_sarar);

比较两个模型的对数似然值和 AIC，选更优的那个作为主回归。审稿人看到你主动做了 SARAR 对照，对模型稳健性的认可度会明显提升。

4.3 动态SDM的估计主流程：集中似然与搜索区间

动态 SDM 的完整设定是：

y_t = τy_{t-1} + ρWy_t + ηWy_{t-1} + X_tβ + WX_tθ + μ_i + ε_t

估计前先做变形。把含 ρWy_t 的项移到左边，得到 (I - ρW)y_t = τy_{t-1} + ηWy_{t-1} + X_tβ + WX_tθ + μ_i + ε_t。这样处理后，右侧的 W 矩阵不再直接作用于 y_t，而是并入设计矩阵：

% 1. 构造滞后变量和空间滞后 ylag = lag_series(y, N, T); % y(t-1)，需要每地区内部滞后 wylag = W * ylag; % W y(t-1) % 2. 设计矩阵 D 不含 ρWy 项 % (I - ρW)y = τ*y_lag + η*W*y_lag + X*β + WX*θ + ε D = [ylag, wylag, x, W * x]; % 3. 计算 rho 搜索区间 lambda = eig(W); rho_min = 1 / min(lambda) + 1e-4; rho_max = 1 / max(lambda) - 1e-4; % 4. 一维搜索 rho：f2_sar_jihai 是集中对数似然函数 opts = optimset('TolX', 1e-10); rho_hat = fminbnd(@(r) -f2_sar_jihai(r, y, D, W), rho_min, rho_max, opts); % 5. 用最优 rho 回代得到全部系数 results = f2_sar_jihai(rho_hat, y, D, W); prt_sardynamic(results);

逻辑说明：f2_sar_jihai 内部先按给定的 ρ 计算 y_star = y - ρWy，然后对 y_star = D·δ + ε 做 OLS，得到残差平方和并计算集中对数似然。外层 fminbnd 只搜索 ρ 一个维度，收敛速度快很多。

参数说明：rho_min 和 rho_max 必须严格落在特征值区间内。上面对 max(lambda) 的统一写法有个前提是 λ_max = 1（行标准化后）。如果 W 已经标准化，rho_max 直接取 0.999 就安全。但我在例行检查里仍然会显式计算特征值，防止 Excel 里混入了未标准化的行列。

动态 SDM 比静态更难收敛，因为多了 τ 和 η，似然面更平坦。常见做法是先跑静态 sar_jihai 拿到 rho0，然后以 rho0 附近作为动态模型的搜索起点：

% 以静态结果为初值，只在小范围内搜索动态版本 static_result = sar_jihai(y, W, x, T); rho0 = static_result.rho; rho_opt = fminbnd(@(r) -f2_sar_jihai(r, y, D, W), ... max(rho_min, rho0 - 0.1), ... min(rho_max, rho0 + 0.1), opts);

这样能避开不少局部最优。等跑完再把 rho_opt 与 rho0 对比：如果两者差到 0.1 以上，说明动态项的加入改变了空间溢出结构，需要回看模型设定是否合理。

5. 避坑指南：动态空间面板回归里最常见的5个翻车现场

5.1 现象：rho输出几乎等于1，似然值异常

最常见的一种翻车是把 rho 估计顶到区间边界，rho 输出接近 0.999 甚至等于 1，对数似然值大得离谱。最初我以为模型设定错了，反复调变量，后来检查发现是权重矩阵没做行标准化。行和不等于 1 意味着 W 的最大特征值远超 1，rho 的可行域整体偏移，优化器直接把 rho 推到了搜索边界上。

解决：每次读入 W 后强制走标准化检查。我现在把这段写成固定防御代码：

if abs(max(sum(W, 2)) - 1) > 1e-10 W = W ./ sum(W, 2); warning('W 已重新做行标准化'); end

从那之后，只要 rho 顶边界，先查 W，再查搜索区间，基本一次定位。

5.2 现象：动态项τ的符号反转，η显著为负且与直觉矛盾

在某个产业数据上，我本来预期时间惯性系数 τ 为正，结果跑出来显著为负；同时 η（时空扩散）也显著为负，整个结果解释不了。问题是 y_{t-1} 和 Wy_{t-1} 高度相关。产业数据里相邻省份的产出往往同涨同跌，滞后值与空间滞后值的相关系数能达到 0.9 以上，一起放进方程后，两个系数的标准误被放大，符号自然不稳定。

解决：先查相关矩阵再决定是否同时保留 τ 和 η。如果相关系数超过 0.85，建议只保留 η 去掉 τ（或者相反），改用动态 SAR 版本：

% 共线性诊断 lag1 = lag_series(y, N, T); corr_matrix = corr([y, lag1, W * lag1, x]); disp(corr_matrix); % 如果 lag1 与 W*lag1 相关系数过高，去掉 η 项 D_sub = [lag1, x, W * x]; % 只保留时间滞后，不加空间滞后

审稿人看的是模型稳健性，不是变量堆得越全越好。两个高度相关的滞后项硬塞在一起，只会让结果看起来更不可信。

5.3 现象：f2_sar_jihai报“Matrix is singular”

报这个错，通常是 W 里有全零行——某个区域没有任何邻居。常见于用经济距离矩阵时，某个特别偏远或经济体量极小的地区被阈值为 0 切断，整行都是 0，空间乘子矩阵求逆时直接奇异。

解决：定位全零行的位置，决定是删样本还是手动补邻居：

% 找出全零行 zero_rows = find(sum(abs(W), 2) == 0); if ~isempty(zero_rows) fprintf('全零行：%s\n', num2str(zero_rows')); % 方案1：手动指定最近邻 for i = 1:length(zero_rows) r = zero_rows(i); [~, j] = min(dist_mat(r, :)); % dist_mat是地区间距离矩阵 W(r, j) = 1; W(j, r) = 1; end W = W ./ sum(W, 2); end

逻辑说明：手动补邻居的前提是你知道这个地区在真实经济网络中应该和谁关联，而不是盲目填 1。如果实在无法确定，删掉该样本是更稳妥的做法。

5.4 现象：xlsread返回的矩阵出现大面积NaN

跑一半发现 y 或 X 读进来全是 NaN，第一反应是数据类型有问题，但实际往往更简单——Excel 文件正被打开，目录下残留了“~$副本高技术产业生产率与集聚W.xlsx”这类临时锁定文件，MATLAB 读取时锁冲突导致部分数据读不出；另一个常见原因是 Excel 单元格里混入了中文表头，默认从第 1 行读取时把文本列读成了 NaN。

解决：先关掉所有打开的 Excel，删除目录下以“~$”开头的临时文件，然后明确指定读取区域：

% 关Excel后，确认没有~$文件残留 dir('~$*'); % 从第2行开始读数据，第一行作为变量名 [data, txt] = xlsread('data.xlsx', 1); % 或指定sheet

读取后用 isfinite 做一次全矩阵体检：任何一个变量列出现 NaN，立刻在数据源排查，而不是在模型里补缺失值。

5.5 现象：longterm.m输出与手算结果对不上

跑完长期效应，longterm.m 输出的数值和我用公式手算的完全对不上。原因在于长期效应要用空间乘子矩阵 (I - ρW)^(-1) 展开，不能拿普通非空间面板的长期乘子公式直接套。很多人在这一步拿 β/(1-τ) 算长期效应，把空间结构忘得一干二净。

解决：显式计算空间乘子矩阵，验证 longterm.m 的内部逻辑：

% 手算总效应：对第k个解释变量 S = inv(eye(N) - rho * W); beta_k = results.beta(k); theta_k = results.theta(k); % 直接效应 direct_k = mean(diag(S) * beta_k); % 间接效应（溢出效应） indirect_k = mean(sum(S, 2) .* theta_k); % total = direct + indirect total_k = direct_k + indirect_k; fprintf('变量%d：直接效应%.4f，间接效应%.4f，总效应%.4f\n', ... k, direct_k, indirect_k, total_k);

参数说明：diag(S) 取乘子矩阵对角线，对应每个地区的自身传导；sum(S, 2) 累计所有邻居的传导路径，是计算间接效应的核心。手算结果与 longterm.m 一致后，才能放心写进论文表格。

6. 复现一份实证结果：从data.xlsx到论文级输出

把整个过程串起来，我通常用一个脚本从上跑到下，中间任何一个环节报错就停下来修：

% ===== 完整复现流程 ===== % 1. 读数据与权重 raw = xlsread('data.xlsx'); y = raw(:, 3); x = raw(:, 4:end); W = xlsread('matrix.xlsx'); W = W ./ sum(W, 2); % 行标准化 % 2. 构造滞后变量 N = 31; T = 13; % 按实际面板维度改 ylag = lag_series(y, N, T); wylag = W * ylag; % 3. 动态SDM设计矩阵 D = [ylag, wylag, x, W * x]; % 4. rho搜索区间 lambda = eig(W); rho_min = 1 / min(lambda) + 1e-4; rho_max = 1 - 1e-4; % 5. 一维搜索并回代 opts = optimset('TolX', 1e-10); rho_hat = fminbnd(@(r) -f2_sar_jihai(r, y, D, W), rho_min, rho_max, opts); results = f2_sar_jihai(rho_hat, y, D, W); % 6. 输出论文表格 prt_sardynamic(results); prt_sp(results);

跑完之后我固定做三件事验证结果能不能写进论文。第一，确认 rho_hat 不在搜索区间的边界上，贴边跑的模型基本不可信。第二，把模型拟合值画出来和实际值对比：plot(y), hold on, plot(results.yhat)，两者相关系数低于 0.9 说明变量选择有问题，回看 X 和 W 定义。第三，用 longterm.m 和手算乘子矩阵对照长期效应，数字对不上就重新排查。

这套流程看起来笨重，但拦下了我过去至少一半的翻车。从那以后，每次拿到新的空间面板数据，我都强制走一遍：标准化权重 → 特征值体检 → 静态模型打底 → 动态模型加滞后 → 长期效应手算核对。顺序换不得，希望帮到你。

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