简介本资源是一套完整、可直接运行的空间双重差分SDID模型MATLAB实现代码面向空间计量经济学研究者、区域政策评估人员及高年级硕博生专为解决含空间溢出效应的准自然实验因果识别问题而设计。包内共65个文件涵盖49个核心.m函数如TW权重生成、DID虚拟变量自动构建、局部/全局效应分解、面板数据堆叠、稳健估计与Hausman检验等、9个.xlsx数据模板含事件/时期虚拟变量示例、多期面板变量结构及7个.mat实证数据与中间结果文件总容量10.62MB模块划分清晰支持从数据预处理、权重内生构造到多模型自动遴选NLSDIDAll的一站式分析。已有1874人学习下载代码经作者亲测可用包含大量注释与演示脚本如零向量分解、个体效应分离等显著降低空间DID方法在Matlab平台的实现门槛特别适合开展区域政策空间溢出效应实证研究的用户快速上手与方法复现。 空间双重差分SDID的Matlab代码网上能找到的版本很多但真正能跑通、能换到自己数据上直接用的其实不多。最近帮几个师弟师妹调空间双重差分代码发现大家卡住的点高度一致权重矩阵构建逻辑不对、极大似然函数迭代不收敛、结果出来符号跟理论预期相反。这篇文章把我亲测可用的SDID全套Matlab代码整理出来从空间权重矩阵构建、面板数据组织、DID变量生成到极大似然估计一步不落。运行环境是Matlab R2021a以上不依赖额外的空间计量工具箱只用基础函数加优化工具箱就能跑通。如果你正在做政策评估、区域经济、公共财政、环境规制这类涉及空间溢出的实证研究这篇文章刚好能帮你省掉至少两三个星期的试错时间。1. 为什么你的DID模型需要给“空间”留个位置1.1 传统DID的隐形假设政策效应不会“串门”传统双重差分模型的基本形式是Y_it α β·(Treat_i × Post_t) γX_it μ_i λ_t ε_it这个设定看着简单其实背后藏着一个很强的假设个体之间互相独立处理组受到的政策影响不会扩散到对照组。换句话说它默认政策效应被关在“笼子”里只影响本地。但现实中的政策很少有这种“绝缘”属性。举个例子某区域设立经济开发区本地GDP确实涨了但周边地区的企业为了享受配套红利可能也会跟着受益甚至一部分产业会从周边迁入开发区导致周边GDP反而下降。这两种情况都意味着对照组的结果变量被处理组的政策间接改变了。传统DID如果忽略这种空间溢出β估计量会同时包含“政策对处理组的直接影响”和“政策对对照组的溢出效应”两者混在一起没法区分。更麻烦的是误差项也可能存在空间相关导致标准误估计偏低显著性检验虚高。1.2 空间溢出会让你的估计偏多少空间溢出对DID估计的污染方向取决于溢出的符号。如果政策产生了正向溢出也就是处理组受益的同时也带动了周边那么传统DID会高估政策效果因为对照组也变好了相对差异反而变小。如果政策是“虹吸效应”处理组抢走了周边的资源那么传统DID会低估政策效果因为对照组的Y被拉低了处理组和对照组的差异被人为放大。实际研究中更隐蔽的情况是空间相关性主要存在于误差项中。比如两个相邻地区共享一些不可观测的冲击它们之间的误差项相关。这种情况下DID的系数估计量虽然仍然可能是一致的但标准误被严重低估通常会让t值偏大得出去显著但实际不稳健的结论。所以要不要做空间双重差分本质上不是“我喜不喜欢空间计量”的问题而是你的研究设计里政策效应是否可能跨越行政边界传导。如果答案是“可能”那就需要给模型加上空间结构。1.3 SDID的模型设定从SAR、SEM到SDM空间计量模型经过几十年发展形成了一个家族。放到DID框架里常见的有三种模型设定形式核心含义SAR-DIDY ρWY Xβ δD ε被解释变量存在空间溢出SEM-DIDY Xβ δD u, u λWu ε误差项存在空间相关SDM-DIDY ρWY Xβ WXθ δD W_D·θ_d ε解释变量和被解释变量都存在空间交互其中SDM-DID空间杜宾DID是最常用的形式因为它同时容纳了两种效应空间滞后项ρWY刻画的是“邻居的结果会影响本地结果”空间滞后解释变量WX刻画的是“邻居的特征会影响本地结果”。放到DID情境里我们通常还会把DID交互项D的空间滞后WD也放进去用它来识别政策对邻近地区的溢出效应。我个人的习惯是优先跑SDM-DID然后通过LR检验或Wald检验看能不能退化成SAR或SEM。这样做的好处是避免从一开始就设定错误后面还可以给审稿人一个完整的模型选择逻辑。2. 跑SDID之前必须做好的三件事2.1 空间权重矩阵模型的“地基”空间权重矩阵W是SDID最核心的输入它决定了“邻居”的定义。W的每个元素w_ij表示地区j对地区i的影响权重。常见的构造方式有邻接矩阵相邻为1不相邻为0距离矩阵距离的倒数或负指数经济距离矩阵经济特征相似度的倒数在Matlab里构建邻接矩阵是最简单的但对于不规则区域比如行政边界复杂的省份最好还是导入ArcGIS或GeoDa生成的边界信息。这里强调一个几乎所有新手都会忽略的问题W必须做行标准化。每一行求和每个元素除以该行和使得每行之和为1。行标准化后的WY可以解释为“邻居Y的加权平均”ρ直接反映邻居均值对本地的影响强度解释起来非常直观。更重要的原因是数学上的行标准化后的W谱半径最大为1只要|ρ|1矩阵(I−ρW)就是非奇异的极大似然估计的雅可比项才有定义。如果W不标准化ρ的参数空间会变得很不规则估计时很容易撞上奇异矩阵。2.2 面板数据如何嵌入空间矩阵空间权重矩阵W是N×N的而你的面板数据是N×T条记录两者怎么匹配是另一个经典坑点。正确的做法是把W扩展成分块对角矩阵W_panel kron(I_T, W)其中I_T是T×T的单位矩阵kron表示Kronecker积。这样扩展后的矩阵是(NT)×(NT)的每个时间截面上的空间关系都由同一个W描述。也就是说W_panel把N×T的向量按“先个体、后时间”的堆叠顺序作用每个时间切片的邻居关系保持一致。这段逻辑如果理解不了建议直接用代码验证一小步随机生成一个N×T矩阵按stack方式变成长向量再用kron矩阵乘一遍看看是否等于每个时间截面上分别做W乘法。我在第一次写这段代码时就是用这种方式验证的确认无误后心里才踏实。2.3 为什么是极大似然估计空间滞后项ρWY是内生的因为它和误差项ε相关直接用OLS估计会得到有偏且不一致的结果。空间计量文献的常用估计方法是极大似然估计ML或广义矩估计GMM。ML的核心是对数似然函数中包含一个雅可比项T·log|I_N − ρW|它来自误差项从ε到Y的变换也是把ρ识别出来的关键。ML估计的流程不复杂给定参数ρ、β、σ²计算残差e (I − ρW)Y − Xβ代入对数似然函数用数值优化算法搜索使对数似然最大的参数值。难点在于雅可比项的计算和数值优化的稳定性这部分我会在代码注释里详细说明。3. 亲测可用的Matlab全套代码3.1 数据导入与模拟数据生成我先给出一段生成模拟数据的代码作用是让整套流程完整跑通。实际研究时你把这段换成readtable读取自己的Excel数据即可。clear; clc; close all; rng(2025); % 固定随机种子保证结果可复现 % 基础参数 N 30; % 个体数比如30个地区 T 10; % 时间期数比如10年 NT N * T; % 面板总样本量 % 如果跑真实数据从这里替换 % data readtable(你的数据.xlsx); % id列、时间列、Y列、X列都需要你根据自己数据调整 % 下面所有变量生成部分直接用data里的列替换即可 % 生成模拟数据的逻辑是设定一组真实的参数β、ρ、θ、δ然后反解出Y。这样跑完估计后你能直接对比估计值和真实值的差距验证代码是否写对了。这是我最推荐的代码自检方式先用模拟数据确认估计能还原参数再上真实数据。3.2 空间权重矩阵构建构建邻接矩阵时我使用了“距离带”的方式模拟地区相邻关系每个地区和前后两个地区相邻。% 构建空间权重矩阵邻接结构 adj (abs((1:N) - (1:N)) 2) - eye(N); % 这行代码生成了一个带状邻接矩阵i和j的索引差不超过2就认为相邻 W adj ./ sum(adj, 2); % 行标准化每一行除以行和 W(isnan(W)) 0; % 防御性处理避免除零这里sum(adj, 2)是每行求和./(行和)实现了行标准化。你如果用的是GeoDa生成的权重矩阵文件导入后也一定记得标准化。真实研究中的邻接矩阵多半不是这种规则结构。比如研究中国地级市的话常见做法是导入地级市边界shp文件在ArcGIS或GeoDa里生成queen邻接矩阵导出成矩阵或csv后读入Matlab。无论哪种方式W的维度必须是N×N且W(i,j)表示j对i的影响别把方向搞反了。3.3 生成DID交互项与空间滞后项DID交互项D的构造逻辑和传统DID完全一致处理组虚拟变量乘以政策后时期虚拟变量。% 生成DID交互项 treat zeros(N, 1); treat(1:15) 1; % 前15个地区为处理组 post zeros(T, 1); post(6:10) 1; % 第6期开始为政策实施后 % 扩展成面板长度treat按个体重复T次post按时间重复N次 D kron(treat, ones(T,1)) .* kron(ones(N,1), post);关键在这里kron(treat, ones(T,1))会生成一个NT×1的向量前15个地区对应的所有时期都是1后15个地区都是0。kron(ones(N,1), post)会让前5期所有地区都是0、后5期所有地区都是1。两者逐元素相乘得到的就是“处理组×政策后”取1的DID交互项。生成空间滞后变量的方法是统一乘W_panel% 面板空间权重矩阵与空间滞后变量 Wbig kron(eye(T), W); % (NT × NT)分块对角矩阵 X1 randn(NT, 1); % 解释变量1模拟数据真实数据直接替换 X2 randn(NT, 1); % 解释变量2 WX1 Wbig * X1; % X1的空间滞后 WX2 Wbig * X2; WD Wbig * D; % DID交互项的空间滞后这是SDID的关键项这段代码里Wbig * X1的含义是对每个时间截面分别计算空间权重加权平均。因为Wbig是分块对角的乘一次等于对所有时期同时做空间加权效率很高。3.4 对数似然函数实现对数似然函数是整段代码的灵魂。我用了tanh变换把ρ限定在(−1,1)区间内避免优化过程中ρ越界导致行列式非正。function nll sdid_nll(paras, Y, X, W, N, T) % paras [atanh(rho); beta; log(sigma)] % 使用tanh变换保证rho在(-1,1)内log变换保证sigma为正 rho tanh(paras(1)); % 还原rho beta paras(2:end-1); % 回归系数向量 sigma exp(paras(end)); % 还原sigma % 计算残差 e (I - rho*W_panel) * Y - X * beta Wbig kron(speye(T), W); I_NT speye(N * T); e (I_NT - rho * Wbig) * Y - X * beta; % 雅可比项 T * log|I_N - rho*W| % 用LU分解计算行列式的对数比直接det()更稳定 IN speye(N); A IN - rho * W; [~, U] lu(A); logabsdet sum(log(abs(diag(U)))); % 对数似然 loglik -(N*T)/2 * log(2*pi) ... - (N*T)/2 * log(sigma^2) ... T * logabsdet ... - (e * e) / (2 * sigma^2); nll -loglik; % fminunc是最小化器所以返回负对数似然 end关于雅可比项有两点需要解释。第一为什么是T乘log|I_N − ρW|而不是log|I_NT − ρW_panel|因为W_panel I_T ⊗ WI_NT − ρW_panel I_T ⊗ (I_N − ρW)而分块对角矩阵的行列式等于各块行列式的乘积所以log行列式等于T × log|I_N − ρW|。第二用LU分解而不是直接det(A)是因为det在高维时可能溢出或下溢LU分解后对U的对角线取绝对值再求和数值稳定性好很多。3.5 主程序估计与结果输出主程序部分负责生成Y、设置初值、调用fminunc优化并输出结果。% 生成被解释变量Y模拟数据用 rho_true 0.5; beta_true [1.0; -0.8]; delta_true 1.2; theta_true [0.3; 0.2; -0.5]; % 对应WX1, WX2, WD sigma_true 0.2; % 注意这里为了演示X_mat的列顺序是 [X1, X2, D, WX1, WX2, WD] X_mat [X1, X2, D, WX1, WX2, WD]; beta_all [beta_true; delta_true; theta_true]; % 解线性方程组 (I - rho*W_panel) * Y X*beta epsilon I_NT speye(NT); epsilon sigma_true * randn(NT, 1); Y (I_NT - rho_true * Wbig) \ (X_mat * beta_all epsilon); % 极大似然估计 % 初始值rho0.3, beta先猜0.5, sigma0.5 paras0 [atanh(0.3); 0.5 * ones(size(beta_all)); log(0.5)]; objfun (p) sdid_nll(p, Y, X_mat, W, N, T); options optimoptions(fminunc, ... Algorithm, quasi-newton, ... Display, iter, ... MaxFunctionEvaluations, 2000, ... MaxIterations, 500); [paras_hat, fval, exitflag] fminunc(objfun, paras0, options); % 还原参数并打印 rho_hat tanh(paras_hat(1)); beta_hat paras_hat(2:end-1); sigma_hat exp(paras_hat(end)); fprintf(\n SDID 估计结果 \n); fprintf(rho_hat %.4f (真实值 %.4f)\n, rho_hat, rho_true); fprintf(beta1_hat %.4f (真实值 %.4f)\n, beta_hat(1), beta_true(1)); fprintf(beta2_hat %.4f (真实值 %.4f)\n, beta_hat(2), beta_true(2)); fprintf(delta_hat %.4f (真实值 %.4f)\n, beta_hat(3), delta_true); fprintf(theta1_hat %.4f (真实值 %.4f)\n, beta_hat(4), theta_true(1)); fprintf(theta2_hat %.4f (真实值 %.4f)\n, beta_hat(5), theta_true(2)); fprintf(thetaD_hat %.4f (真实值 %.4f)\n, beta_hat(6), theta_true(3)); fprintf(sigma_hat %.4f\n, sigma_hat); fprintf(Log-likelihood %.4f\n, -fval); fprintf(Exit flag %d\n, exitflag);模拟数据的妙处在于跑完可以直接把估计值和真实值做对照。如果代码正确估计值应该在真实值附近小幅波动。真实数据没有“真实值”可对照你就要靠系数符号、显著性和经济含义来判断。这里要特别说明上面为了聚焦SDID核心机制我没有把个体固定效应和时间固定效应放进模拟过程所以估计能还原参数。真实数据里个体/时间固定效应几乎一定存在必须控制具体做法见第6章。4. 结果怎么解读系数、显著性与效应分解4.1 核心系数的经济含义在SDM-DID设定下有几个系数需要重点关注δ对应D的系数含义是“在控制了空间溢出后政策对处理组本地的直接影响”。这是你论文里最想报告的那个数。θ_d对应WD的系数含义是“邻居是否为处理组、是否处于政策期对本地Y的影响”。如果θ_d显著为正说明政策有正向空间溢出邻居受政策影响的同时拉动了本地。ρ对应WY的系数含义是“邻居Y的加权平均对本地Y的影响”。ρ显著说明结果变量本身存在空间传染性比如GDP增长会带动周边增长环境污染会扩散到邻市。4.2 有空间滞后就不能只看回归系数直接效应与间接效应空间杜宾模型一个让新手容易翻车的地方是β和θ不能直接解释为边际效应。因为模型里有ρWY项(I − ρW)会作用于所有解释变量任何一个变量变化都会通过空间反馈机制扩散到所有地区。假设模型是Y ρWY Xβ WXθ ε整理后Y (I − ρW)⁻¹ (Xβ WXθ ε)所以X_k对Y的偏效应是∂Y/∂X_k (I − ρW)⁻¹ (Iβ_k Wθ_k)这是一个N×N矩阵。对角线元素的均值叫直接效应非对角线元素的均值叫间接效应也就是溢出效应两者之和叫总效应。这才是你应该报告的数字。如果只报β审稿人问一句“你的空间反馈效应呢”你很难答好。4.3 效应分解的Matlab实现基于上节公式效应分解的代码很简单% p a hrefhttps://download.csdn.net/download/weixin_61675495/79337417 stylecolor:#ec7500;font-size:14px; 本文还有配套的精品资源点击获取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p