简介：面向空间计量经济学研究者的动态面板空间杜宾模型Matlab实现资源，适合硕博研究生、科研人员及空间计量课程学习者，适用于需同时处理空间依赖与时间动态效应的实证场景，弥补传统静态空间模型未纳入滞后项的局限。压缩包共19个文件，以.m脚本为主体，辅以xlsx数据文件与PDF说明文档，资源仅337KB，轻量易用。脚本覆盖模型估计、空间权重矩阵构建、示例驱动、结果打印等关键环节，PDF系统梳理了动态空间面板极大似然估计的原理与实现思路，Excel数据可直接用于复现或替换为自有数据。目前已有1949人学习下载。借助这份代码，读者可快速跑通动态空间杜宾模型，理解空间动态回归的设定、估计与输出流程，也能参考脚本结构修改模型变量，支撑论文实证、课题研究或课程作业。

1. 动态面板空间杜宾模型：这份压缩包里装了什么，能帮你解决什么问题

动态面板空间杜宾模型，是空间计量里“难做但实证常用”的那一类。它和静态空间杜宾模型的核心区别在于：模型同时允许被解释变量受上一期自身、本期邻居、以及上一期邻居的共同影响，多出来的是两条时间维度的传播路径。换句话说，变量 x 的变化不仅影响本地区当期的 y，还会通过空间溢出和时间惯性，影响邻居下一期的 y。这套压缩包提供的正是一份能直接运行的 MATLAB 实现，配套一个“高技术产业生产率与集聚”的实证数据文件，以及一份 MLE 推导笔记。适合两类人：一类是论文研究了静态空间面板后想补动态效应，另一类是刚入门空间回归，想找一份完整代码当脚手架，而不是从零搭模型推导。

2. 模型内核：动态 SDM 的方程、MLE 估计路线与文件分工

2.1 从静态 SDM 到动态 SDPD：多出来的 τ 和 ψ 不是装饰

先看标准静态空间杜宾模型（SDM）的方程，它写的是当期截面关系：

y_t = ρ W y_t + X_t β + W X_t θ + μ + ε_t

这里面只有两个空间交互来源：ρ 捕捉 y 的当期空间溢出，θ 捕捉解释变量 W X 的空间溢出。所有影响都是“瞬时”的，上一期发生了什么，模型并不关心。

动态空间杜宾（SDPD）在右侧多出两项，变成：

y_t = τ y_{t-1} + ρ W y_t + ψ W y_{t-1} + X_t β + W X_t θ + μ + ε_t

多出来的 τ 是纯时间惯性系数，含义是本地上一期的 y 对本地当期 y 的影响，类似动态面板里的自回归项。而 ψ 是时空滞后系数，也叫 spatio-temporal lag，它刻画的是“邻居上一期的 y 如何影响本地当期 y”。这两个系数是动态 SDM 区别于动态 SAR 或静态 SDM 的关键。

为什么说 ψ 不是装饰？因为在静态模型里，邻居的影响只能通过当期 W y_t 进入，影响的持续性被压缩进 ρ 一个参数里。而在真实数据里，一个地区的产出或技术扩散往往要滞后一期甚至多期才传导到邻居，如果不放 ψ，这种滞后溢出会被错误归因成同期溢出，导致 ρ 高估。这就是动态模型修正静态模型最重要的价值。

2.2 为什么这个模型必须走 MLE：空间内生性与动态面板的偏误风险

遇到这种设定，第一反应是想用 OLS 或 2SLS，但实际行不通。OLS 失效有三个直接原因。第一，W y_t 是内生的，因为把方程重新整理成简化形式后，y_t 可以展开成包含所有地区 y 的多期滞后，右边的 W y_t 和当期误差 ε_t 必然相关。第二，y_{t-1} 在固定效应面板里会带来 Nickell 偏误，空间模型把这个偏误进一步放大了。第三，W y_{t-1} 通常和地区固定效应相关，常规差分或组内变换也未必能彻底清理干净。

走极大似然估计（MLE）是空间计量里比较稳妥的做法。核心思路是把模型写成简化形式：

y_t = (I_N - ρW)^{-1}((τ I_N + ψW) y_{t-1} + X_t β + W X_t θ + μ + ε_t)

变换之后 ε_t 近似服从独立正态分布，于是可以写出完整的对数似然函数。这里有个常用的降维技巧：给定空间参数 ρ、τ、ψ 的条件下，β 和 σ² 可以用类 GLS 公式解析求出，所以优化器只需对少数几个空间参数做搜索。这就是 profile likelihood（剖面似然）思路，也是代码里f2_sar_jihai.m这个目标函数干的事：它接收空间参数组合，返回负对数似然值，sar_jihai.m负责装配数据、调优化器、循环迭代。

这套代码与静态 SDM 的 MLE 实现有个明显区别：静态模型只需要搜索 ρ 一个参数，而动态模型至少要搜索 τ、ρ、ψ 三个，搜索维度和似然面复杂度都上去了。迭代初值怎么给，后文实战部分会专门讲。

2.3 压缩包文件怎么分工：一份能自洽的代码骨架

我解压之后把文件列表过了一遍，这套文件不是散件，而是按“估计核心 + 变体 + 工具 + 示例”组织的完整骨架。

这里给一个建议：拿到压缩包先不要急着跑sar_jihai_example.m，先花半小时把note-sdpd-mle.pdf翻一遍。它解释了似然函数每一项的来历，尤其是雅可比项和初始条件处理。后面排错时，90% 的问题都要回到这个笔记里找依据。SARAR.m和longterm.m是副产物，前者适合做对照模型，后者正好是实证论文里最需要的长期效应。

3. 把示例跑通：从数据准备到 sar_jihai_example.m 输出解读

3.1 解压之后先删掉 ~$ 开头的临时文件

我解压这份 rar 之后，第一眼就发现了两个隐藏坑：~$副本高技术产业生产率与集聚W.xlsx和~$副本高技术产业生产率与集聚data.xlsx。这类~$前缀文件是 Office 打开 Excel 时自动生成的临时锁文件，不是有效数据。用 WinRAR 或 7-Zip 解压后，如果不清理直接让 MATLAB 去读，xlsread在 Windows 下可能报“文件正由另一进程使用”，或者读入脏数据。

我的习惯是先做一次清理再动数据：

cd /path/to/unzipped rm -f '~$*' ls -la

清理后目录里剩下的就是真正的数据文件。另外注意检查 rar 包是否设置了只读属性，解压后先把属性改成可写，否则 MATLAB 写临时文件时会失败。这一步虽然简单，但能省掉后面不少莫名其妙的报错。

3.2 入口脚本的典型结构：读数据、调权重、估计、输出

sar_jihai_example.m是整个包的入口。它做的事情可以拆成四段，每一段对应一个明确目标。下面按这类 MATLAB 代码最常见的结构还原出来，便于对照着看：

%% 1. 读入数据 raw = xlsread('data.xlsx'); % 面板数据：id, year, x1, x2, ..., y id = raw(:,1); % 地区编号 year = raw(:,2); % 年份 X = raw(:, 3:end-1); % 解释变量矩阵 y = raw(:, end); % 被解释变量 %% 2. 读权重矩阵并行标准化 Wc = xlsread('matrix.xlsx'); % N x N 邻接矩阵 W = wrook('matrix.xlsx'); % 行标准化，得到 W %% 3. 设置初值并估计 T_data = length(unique(year)); N = size(W, 1); theta0 = [0.5, 0.2, 0.1]; % 初值顺序 tau, rho, psi options = optimset('Display','iter', 'MaxFunEvals', 5000, 'TolX', 1e-8); [theta_hat, nll] = sar_jihai(y, X, W, T_data, theta0, options); %% 4. 打印结果 prt_sardynamic(theta_hat, X, W, y, nll);

这里参数顺序很重要：theta0我按[tau, rho, psi]三个空间参数的顺序给。初值不要随便给 0，常见做法是先跑一遍不带动态项的静态 SDM，把它的 ρ 估计值作为rho的初值，tau取 0.5 左右，psi取 0 或rho/2。这样优化器起步位置离真实值更近。

TolX设到 1e-8 也是有意为之。动态空间模型的似然面往往比静态模型平坦很多，默认的 1e-6 可能停在鞍点上，迭代日志里看到参数不再变化但似然值还在小幅波动时，先收紧这个容差再重跑。

3.3 输出结果表怎么读：哪些系数直接解释，哪些不能

prt_sardynamic.m输出的表格里，核心看这几组内容：

第一，三个关键系数：rho显著为正代表同期限空间溢出存在，这是静态模型里也有的结论；tau显著为正说明存在时间惯性，地区当期表现受上期自身影响；psi显著则说明确实存在“邻居上期 -> 本地当期”的滞后溢出路径，这是动态模型的核心证据。第二，解释变量系数 β 和它们的空间滞后项系数 θ。第三，对数似然值、R² 之类拟合指标。

这里要专门提醒一句：动态模型里的 β 不能直接当边际效应解释，因为 W y_t 和 W y_{t-1} 的存在会让一个解释变量变化引发的冲击在空间和时间两个维度上反复传播。真正能写在论文表格里的边际效应，是短期直接效应、短期间接效应和长期效应，下一章专门讲怎么算。如果你看到别人的论文只报 β 不报效应分解，基本可以判断他没做完整。

4. 动态效应分解：短期与长期乘子，才是动态 SDM 的真正输出

4.1 效应分解的数学准备：直接效应与间接效应从哪来

空间计量里谈边际效应，用的是 LeSage 和 Pace 定义的偏导数矩阵方法。对某个解释变量 k，先构造矩阵 S_k = (I_N - ρW)^{-1}(β_k I_N + θ_k W)，然后取它的对角线平均值作为直接效应，取行和减去对角线平均作为间接效应。直接效应衡量本地区解释变量变化对本地区 y 的影响，间接效应衡量它对所有其他地区 y 的影响，两者相加是总效应。

动态模型比静态多一层麻烦：要区分“短期”和“长期”。短期看的是当期冲击对当期 y 的影响，直接套上面的 S_k 即可；长期则要考虑时间上的累积传播，把 y 的通解演变到稳态。

在稳态下 y_t = y_{t-1} = Y*，把动态方程代进去整理，可以得到：

((1-τ)I_N - (ρ+ψ)W) Y* = Xβ + W X θ

所以长期乘子矩阵是：

M_k = ((1-τ)I_N - (ρ+ψ)W)^{-1}(β_k I_N + θ_k W)

注意这里矩阵求逆的相对位置和 S_k 不同，区别来自时间维度的累积。这也是 dynamic SDM 和静态 SDM 效应分解最本质的区别：如果误用 S_k 当长期效应，会把同时期溢出与动态累积混为一谈。

4.2 用 MATLAB 还原效应计算脚本

把上面的公式翻译成 MATLAB 代码，逻辑并不复杂，关键在于别把矩阵运算方向搞错。下面是一段可以直接套用的效应计算脚本：

function [di_s, si_s, di_l, si_l] = sdpd_effects(beta, theta, rho, tau, psi, W, N) % beta 和 theta 是某个解释变量的系数组合 % 返回短/长期直接效应与间接效应 A_s = eye(N) - rho * W; % 同期限空间变换 S_k = A_s \ (beta * eye(N) + theta * W); % 短期乘子矩阵 di_s = trace(S_k) / N; % 短期直接效应 si_s = (ones(1,N) * S_k * ones(N,1)) / N - di_s; % 短期间接效应 A_l = (1 - tau) * eye(N) - (rho + psi) * W; % 长期均衡矩阵 M_k = A_l \ (beta * eye(N) + theta * W); % 长期乘子矩阵 di_l = trace(M_k) / N; si_l = (ones(1,N) * M_k * ones(N,1)) / N - di_l; end

代码说明：trace(S_k)是 LeSage 直接效应的标准公式；ones(1,N) * S_k * ones(N,1)计算的是矩阵所有元素之和，等价于总效应的 N 倍，减去对角线就是全部间接效应。A_s \是 MATLAB 的左除运算，求解线性方程组，比写inv(A_s) * ...数值上更稳定，尤其当 N 到了几百个地区时。

有一个需要盯住的边界：当(rho + psi) / (1 - tau)趋近 1 时，长期矩阵A_l接近奇异，效应值会爆炸式增长。这说明模型认为地区间的相互依赖已经接近不稳定状态，这时候要先检查估计结果，再考虑是否模型设定有问题。跑完效应分解，最好顺手把cond(A_l)打出来看一眼。

5. 避坑手册：空间权重矩阵和数据排序引发的五个翻车现场

5.1 权重行列顺序与面板排序不一致

现象：估计结果里对数似然值极差，甚至出现负的 R²，ρ 的估计值和用静态模型得到的参考值完全对不上。

原因：data.xlsx里的数据行顺序和matrix.xlsx的行顺序不是同一个排序。权重矩阵 W 第 i 行对应地区 i，面板数据如果按年份堆叠，行号不一定和 W 的行号对应。一旦错位，W * y乘出来的就是“张冠李戴”的空间滞后项。

解决：读数据后先排序，再调权重。我每次都会写一个自检片段：

unique_id = unique(id); if length(unique_id) ~= size(W, 1) error('面板地区数与权重矩阵维度不一致'); end % 核对 W 的行顺序与 unique_id 完全一致后再进入估计

这个检查必须放在wrook调用之前，因为行标准化不会修正顺序错误，只会把错误顺序也标准化得“看起来很合理”。

5.2 W 矩阵包含自环或孤立点

现象：某个地区的所有估计效应都是 0，或者lambda_max特征根算出来大于 1，模型不稳定。

原因：邻接矩阵对角线上出现了 1，也就是某个地区把自己当邻居；或者某一行全为 0，规范标准化后该行仍然全 0，W * y对那个地区恒为 0。

解决：读入 W 后先做诊断。常见做法是检查diag(W)是否全 0、sum(W,2)是否有 0 行。对于孤立点，要么剔除该地区，要么在文献里说明并按最小邻接处理，不能放任不管。代码里加一句：

if any(sum(W,2) == 0) warning('存在孤立点地区，请检查邻接定义'); end

5.3 动态项和时空滞后项高度共线

现象：τ 和 ψ 的符号相反、标准误非常大，似然面像一条长长的沟谷，迭代半天参数不动。

原因：y_{t-1}和W y_{t-1}高度相关。如果某个地区的 y 在空间上有强集聚性，上一期自身的值和上一期邻居的值几乎线性重合，这时 τ 和 ψ 很难被分别识别。

解决：先从只含 ρ 和 τ 的简化模型入手，确认动态项确实有增量价值，再加 ψ。我一般会用 BIC 做一次比较，而不是直接上完整模型。如果加入 ψ 后 logL 只提高一点点而 τ、ψ 标准误暴涨，基本可以判断共线，老老实实报告简化版本。

5.4 时间维度 T 太小导致似然面平坦

现象：迭代不收敛，或者 ρ 的 t 统计量低到不合理；无论初值怎么换，结果都在某个范围内漂移。

原因：动态面板的初始条件问题在空间模型里一点没少。当 T 只有 5～10 期而 N 很大时，MLE 的偏误明显，似然函数对 τ、ψ 的识别力不足，这是模型的数学性质决定的，不是代码 bug。

解决：T 小于 10 的数据，优先考虑用 Lee 和 Yu 那类偏误修正估计量；或者报告结果时明确标注“动态项为探索性分析”。拿到这套代码后，先确认note-sdpd-mle.pdf里的推导是否包含了偏误修正项，如果没有，那它更适合长面板场景。

5.5 rar 解压陷阱与临时文件混入

现象：xlsread读出的变量数量比预期多，或者报错“文件正在使用”。

原因：压缩包里混入了~$开头的 Office 临时锁文件，还有以~$副本命名的文件。MATLAB 在做dir('*.xlsx')时会把它们一并匹配进来，在 Windows 下就会触发读取冲突。

解决：解压后第一件事就是删除所有~$*文件。这一步在 3.1 节里已经提过，但值得再强调一次，因为这是最容易在项目初期浪费半小时的玄学问题。另外，如果 rar 包是从网盘下载的，先检查完整性再解压，避免读到一半报 CRC 错误。

6. 把代码改造成自己的研究：数据替换、参数调整与可信度自查

6.1 新面板数据预处理三步

换到自己数据集时，先做三步标准化操作。第一步，用sortrows按地区和年份排序，保证同一地区的所有年份相邻；第二步，把权重矩阵的行顺序校准成排序后的地区顺序；第三步，做一次平衡面板检查，看是不是每个地区都有完整的 T 期数据：

[~, idx] = sortrows([id, year], [1, 2]); y = y(idx); X = X(idx, :); tabulate(id); % 检查是否平衡面板

6.2 修改控制变量与调用参数

data.xlsx换成自己的数据后，X 矩阵的列数会变，其余几乎不用动。如果你引入时间固定效应，调用处把sar_jihai换成sar_jihai_time，两者入参格式一致。改完先跑一个简化版本，确认对数似然能正常下降，再逐步加变量。参数初值、优化容差沿用第 3 章的设置即可。

6.3 与静态 SDM 对照：动态项的增量价值与稳健性判断

我的习惯是把静态 SDM 和动态 SDPD 两个结果并排放，比较 logL 和三个核心系数的变化。如果动态版本 logL 显著上升，且 τ、ψ 至少有一个显著，说明动态结构确实存在；如果 logL 只提高了一丁点，τ、ψ 都不显著，那动态项就是自我安慰，不如直接报告静态结果，更诚实也更稳。从那以后，我每次换新数据都强制自己先跑一遍“排序自检 → 权重诊断 → 静态对照”，再谈任何动态结论。空间面板模型的坑大多来自数据整理而不是模型本身，把这个流程走完，代码包才能真正变成你自己的工具。希望帮到你。

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