简介面向空间计量经济学研究者的动态面板空间杜宾模型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.m动态空间杜宾模型主体估计函数f2_sar_jihai.m负对数似然目标函数供优化器迭代调用sar_jihai_time.m加入时间效应/时间趋势的变体f2_sar_jihai_time.m时间变体对应的目标函数sar_jihai_example.m一键运行的示例脚本SARAR.m空间自回归 空间自相关误差的对照模型wrook.m读取邻接矩阵并做行标准化longterm.m长期乘子效应计算prt_sardynamic.m / prt_sp.m / prt_reg.m结果格式化输出note-sdpd-mle.pdf模型推导笔记建议最先读data.xlsx / matrix.xlsx实证数据与空间权重矩阵这里给一个建议拿到压缩包先不要急着跑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规范标准化后该行仍然全 0W * y对那个地区恒为 0。解决读入 W 后先做诊断。常见做法是检查diag(W)是否全 0、sum(W,2)是否有 0 行。对于孤立点要么剔除该地区要么在文献里说明并按最小邻接处理不能放任不管。代码里加一句if any(sum(W,2) 0) warning(存在孤立点地区请检查邻接定义); end5.3 动态项和时空滞后项高度共线现象τ 和 ψ 的符号相反、标准误非常大似然面像一条长长的沟谷迭代半天参数不动。原因y_{t-1}和W y_{t-1}高度相关。如果某个地区的 y 在空间上有强集聚性上一期自身的值和上一期邻居的值几乎线性重合这时 τ 和 ψ 很难被分别识别。解决先从只含 ρ 和 τ 的简化模型入手确认动态项确实有增量价值再加 ψ。我一般会用 BIC 做一次比较而不是直接上完整模型。如果加入 ψ 后 logL 只提高一点点而 τ、ψ 标准误暴涨基本可以判断共线老老实实报告简化版本。5.4 时间维度 T 太小导致似然面平坦现象迭代不收敛或者 ρ 的 t 统计量低到不合理无论初值怎么换结果都在某个范围内漂移。原因动态面板的初始条件问题在空间模型里一点没少。当 T 只有 510 期而 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 只提高了一丁点τ、ψ 都不显著那动态项就是自我安慰不如直接报告静态结果更诚实也更稳。从那以后我每次换新数据都强制自己先跑一遍“排序自检 → 权重诊断 → 静态对照”再谈任何动态结论。空间面板模型的坑大多来自数据整理而不是模型本身把这个流程走完代码包才能真正变成你自己的工具。希望帮到你。本文还有配套的精品资源点击获取