简介这套代码为空间双重差分模型提供了完整的MATLAB实现适合空间计量经济学研究者、研究生及需要开展政策评估实证分析的学者使用。代码围绕模型建模全流程组织涵盖内生时空权重矩阵生成、事件虚拟变量与时期虚拟变量构造、被解释变量与解释变量的面板堆积以及模型估计后的直接效应与间接效应分解并给出基于不同初始权重矩阵的遴选方案便于使用者掌握空间双重差分的关键步骤和参数设定。压缩包共包含六十五个文件其中四十九个脚本、九个表格、七个矩阵数据文件压缩后大小为十点六二兆。脚本按功能模块命名清晰配合示例数据和中间变量文件可在软件中直接运行对照学习。目前已有一千八百七十八人学习使用对于希望快速上手的研究人员能够省去自行整理代码和数据的繁琐过程直接从示例出发修改应用是一条高效的学习路径。1. 空间双重差分SDID的 matlab 代码难的不是语法而是模型里多了两个空间项空间双重差分SDID的 matlab 代码和普通 DID 最大的差别不在写代码的姿势而在你要多处理两样东西空间权重矩阵 W以及由它派生出来的内生空间滞后项 ρWy 和政策交互项 θWD。很多跑 SDID 翻车的人并不是程序报错而是把 W 当成可有可无的装饰结果 ρ 和 θ 全被固定效应吸收D 的系数解释不了政策溢出。这篇笔记按我自己的落地顺序写先讲清楚模型里每一项在识别什么再给出从面板数据到空间权重矩阵、到集中似然 MLE 的可直接改写的 matlab 代码最后是亲测踩过的五个坑。适合已经会用 DID 做政策评估、现在想把溢出效应拆出来的人。读完你可以照框架改自己的数据也能判断拿到的 SDID 结果到底能不能信。2. 模型设定先钉死SDID 每一项到底在识别什么普通 DID 的核心假设是 SUTVA即一个个体的潜在结果不受其他个体处理状态的影响。放到区域政策评估里这个假设几乎总是站不住一个城市搞了产业转移补贴隔壁城市承接迁出企业一个省份开征了环境税下游省份的水质也跟着变。处理效应会穿过地区边界继续传导这就是空间溢出。如果忽略溢出DID 估计出来的系数既不是处理地区的平均处理效应也不是纯对照组的反事实而是两者的混合体。空间双重差分要做的事情就是把这个混合体按来源拆开哪些是处理地区自身的效应哪些是“邻居被处理”传导过来的效应。典型的 SDID 设定在普通 DID 方程上加了两个空间项写成y_it ρ ∑_j w_ij y_jt β D_it θ ∑_j w_ij D_jt X_it γ μ_i λ_t ε_it其中 ρWy 是内生空间滞后意味着本地结果会被同期邻居结果影响这是需要靠极大似然或空间两阶段估计解决的问题WD 是政策交叉项它的系数 θ 是本方法最关心的信号当我周边地区开始被处理时我的结果变量平均改变多少。这两个参数一配合才能在后面的矩阵运算里把总效应拆成直接效应和间接效应。别小看这一步很多人的 SDID 跑出来 D 系数和普通 DID 一模一样就是因为 θ 项根本没进模型固定效应把外溢效应吃掉了。2.1 无干扰假设什么时候站不住对照组被污染的现场先看一个最常见的场景某省推行开发区试点试点城市 A 得到政策邻近城市 B 没有被试点但 A 的招商引资优惠会把原本要落到 B 的企业吸走。此时 B 虽然名义上是对照组它的实际产出路径已经被 A 的处理状态改变了。普通 DID 里 B 的“反事实”不再是反事实而是“邻居被处理后世界下的 B”。于是处理效应估计值偏大还是偏小取决于溢出方向是替代还是互补你根本不知道 β 里混了多少溢出成分。要让 SDID 能识别光加一项 WD 还不够。识别策略的核心在于 W 的结构必须存在一批“既没被处理、邻居也没被处理”的纯对照个体它们为模型不受污染的基准同时还要有一批“自己没被处理、但邻居被处理”的个体它们用来识别 θ。如果政策一旦铺开就几乎覆盖全域W 矩阵又用的是全连接反距离权重那么几乎每个观测都有被处理的邻居参数识别就退化成依靠函数形式假设。这也是为什么我一般在建模前先数一数每种邻居模式的个体数量样本里纯对照太少后面的行列式和标准误再漂亮也没人信。2.2 三种常见空间结构与 SDID 搭配SAR、SEM、SDM 怎么选先把三种常见的空间计量结构说清楚。SAR空间滞后模型只放 ρWy认为结果之间的互相影响是主导机制SEM空间误差模型把空间相关放在扰动项里本质上是处理遗漏空间相关变量的问题SDM空间杜宾模型同时放 ρWy 和 WX即自变量的空间滞后。SDID 的政策设定——也就是 WD 项——其实就属于 WX 的一个特例所以最自然的搭配是 SDM 框架。模型结构方程要点政策识别的含义常见适用场景SARy ρWy Xβ ε结果变量相互影响房价、空气质量这类存在直接交互的变量SEMy Xβ u, u λWu ε空间相关来自共同冲击遗漏重要空间变量时的稳健性检验SDMy ρWy Xβ WXθ ε溢出既通过结果也通过解释变量传导政策评估里同时关心处理和外溢的 SDID 设定实际落地时我一般先跑一个包含 D、WD 以及常见 X 的 SDM 全模型再用似然比检验看 WX 那一块能不能删。能删就退到 SAR 风格简化模型不能删就说明控制变量的空间溢出也存在强行删会让 ρ 偏高。这里有一个容易踩的认知误区加 WD 不是把“政策虚拟变量”简单换成“空间平均的政策虚拟变量”就够而是要在同一套似然函数里同时估计 ρ 和 θ因为 WD 中的 W 是行标准化矩阵WD 会通过 ρ 间接再次影响 y。换句话说直接效应和间接效应都不等于某个单个系数必须按第 4 章的方式展开成 (I - ρW)⁻¹ 的函数。3. 数据准备从面板读到 W 矩阵的完整前奏SDID 对数据格式的敏感度比普通 DID 高很多因为空间权重矩阵要和你面板的堆叠顺序严格对齐。这一章先把数据排列、W 矩阵构造、标准化这三步理顺。代码都是用 matlab 写的常见做法是全部放在一个工作目录里数据文件统一用英文列名避免后面编码问题。3.1 先把面板排成“时间外层、个体内层”的长表用 matlab 处理面板推荐直接用 readtable 读 CSV 或 Excel然后按两列排序。排序顺序直接决定后面 kron 函数的块结构方向这也是初学者最容易翻车的地方。% 读取面板数据CSV 列名: id, time, y, D, X1, X2 data readtable(policy_panel.csv, PreserveVariableNames, true); % 关键一步外层按 time 排内层按 id 排 data sortrows(data, {time, id}); % 把 id 转为数值型城市名这类字符串列必须先 grp2idx if ~isnumeric(data.id) [~, ~, data.id] grp2idx(data.id); end id double(data.id); time double(data.time); y data.y; D data.D; X [data.X1, data.X2]; % 控制变量矩阵 N length(unique(id)); T length(unique(time)); NT N * T;逻辑说明sortrows(data, {time,id})的意思是先把 time 作为主排序键再在同一个时期内部把 id 排好。最终向量里的顺序是 t1 期的 N 个观测接着 t2 期的 N 个观测这种次序称为“时间外层、个体内层”。后面构造 Wbig kron(eye(T), W) 时就是按这个次序生成块对角矩阵如果你这里排的是个体外层那 kron 的参数就要反过来。参数说明grp2idx会把字符串城市名转成 1 到 N 的编号返回的第三个输出可以直接覆盖 data.id这样 id 列就变成 double。T 和 N 的数值要确认和真实面板一致如果数据是平衡面板还好非平衡面板下面临的处理要复杂得多我会先补平衡再用本代码。3.2 用经纬度构造距离权重矩阵K 近邻反距离法地理权重矩阵里最常见的三类K 近邻权重、反距离权重、Queen 邻接权重。如果你手上只有城市的经纬度坐标可以直接在 matlab 里自己造矩阵。function W w_knn_weight(lat, lon, K) % 输入纬度列向量、经度列向量、近邻个数 K % 输出行标准化的 K 近邻反距离权重矩阵 N length(lat); % 经纬度直接算欧氏距离会有误差建议先投影成平面坐标 % 这里用近似平面坐标代替精确做法是 projfwd 投影 xy [lon(:), lat(:)]; D pdist2(xy, xy, euclidean); % 排除自身取最近 K 个邻居的索引 [Dk, idx] mink(D, K1, 2); % 多取一个第 1 个是自身 W zeros(N); for i 1:N ne idx(i, 2:K1); % 去掉自身 dist_ne Dk(i, 2:K1); dist_ne(dist_ne 1e-6) 1e-6; % 避免除零 W(i, ne) 1 ./ dist_ne; % 反距离权重 end % 行标准化放在主函数处理本函数只返回未标准化矩阵 end逻辑说明pdist2计算两两欧氏距离mink(D, K1, 2)表示对每一行取最小的 K1 个值返回距离矩阵 Dk 和对应列索引 idx。由于每行第一列必然是自身所以从第 2 列开始截取就得到 K 个邻居。权重取距离倒数距离越近影响越大。参数说明K 的经验范围是 4 到 10K 太小会让 W 过于稀疏识别出的溢出效应噪声大K 太大则容易让处理组和对照组混在一起θ 变得不显著。另外如果城市分布不均匀比如东部密集西部稀疏固定 K 的权重会让西部城市连接到非常远的邻居这时候可以考虑改用阈值半径法即距离小于某个阈值的地区才相连。3.3 行标准化和孤立点的两步处理这一步省不得空间权重矩阵拿回来之后第一步是行标准化第二步是处理孤立点。行标准化的作用是让 W 的每一行和为 1这样 W 乘以一个变量得到的是“邻居变量的加权平均”ρ 的取值也可以解释为空间溢出强度。% W 来自上文 w_knn_weight 或其它来源 W sparse(W); % 先转稀疏节省内存 % 第一步记录零行 rowsum sum(W, 2); zero_rows find(rowsum 0); % 第二步处理孤立点常见的做法是把自己设为唯一邻居 if ~isempty(zero_rows) for i zero_rows W(i, i) 1; end end % 第三步行标准化 rowsum sum(W, 2); W W ./ rowsum; % 检查是否还有 NaN assert(~any(isnan(W(:))), W 存在 NaN请检查零行处理);逻辑说明零行代表该地区在 K 近邻或半径阈值下没有邻居。如果保留零行行标准化会出现 0/0 的 NaN后面的特征值计算和极大似然全部会报错。把孤立点设成自环的意义是让该地区在空间意义上“只受自己影响”代价是它的空间滞后项退化为自身对 ρ 的识别贡献变小。参数说明sparse(W)在 N 超过 500 时收益明显kron 之后是整个 NT×NT 的大矩阵不用稀疏存储很容易撑爆内存。断言语句是保险丝一旦在调试中触发优先回查零行而不是把矩阵元素改成随机值。这一步看起来简单实际项目里大部分空间权重矩阵的问题都出在“没检查零行”上。4. 核心估计中心化 集中似然 MLE 一次跑通这一步是整份代码的心脏。SDID 里既有面板固定效应又有内生空间滞后项不能像普通 DID 那样直接回归常见做法是用两步 demean 把 μ_i 和 λ_t 消掉然后对剩余方程做极大似然估计。这个方案的优点是不需要生成 NT 个虚拟变量省内存且收敛快缺点是对中心化顺序很敏感顺序错了整个估计结果都是虚的。4.1 两步 demean哪些变量要做顺序为什么是“先中心化后乘 W”对面板数据做个体-时间双向中心化matlab 里用 accumarray 写起来很快。重点在于顺序先把 y、D、X 分别中心化再拿中心化后的变量去乘 W而不是先把 W 乘上去再中心化。因为 demean 矩阵和 W 不可交换先乘 W 会把空间滞后项里混入固定效应的残余成分ρ 容易被高估。function [ydm, Xdm, Ddm] demean_panel(y, X, D, id, time) % 双向固定效应 demean个体均值 时间均值再加回总均值 n length(y); % 对 y 做中心化 id_mean accumarray(id, y, [], mean); time_mean accumarray(time, y, [], mean); grand mean(y); ydm y - id_mean(id) - time_mean(time) grand; % 对 D 做中心化 id_mean accumarray(id, D, [], mean); time_mean accumarray(time, D, [], mean); grand mean(D); Ddm D - id_mean(id) - time_mean(time) grand; % 对 X 每一列做中心化 Xdm zeros(size(X)); for k 1:size(X, 2) id_mean accumarray(id, X(:, k), [], mean); time_mean accumarray(time, X(:, k), [], mean); grand mean(X(:, k)); Xdm(:, k) X(:, k) - id_mean(id) - time_mean(time) grand; end end逻辑说明accumarray(id, y, [], mean)的作用是把相同 id 的所有观测分成一组求组内均值返回一个 N×1 向量再用id_mean(id)把这个均值广播回每个观测。对时间维度同理。最后加回 grand 是因为同时减去个体和时间均值会把总均值减掉两次加回一次保持恒等式完整。参数说明这个函数只适用于平衡面板非平衡面板的 accumarray 分组逻辑没有变化但每个个体贡献的时期数不同均值中心化后数据量仍保持一致只是解释上更复杂。注意 Xdm 的循环是对列进行的如果 X 有几十列循环不会慢因为每列都是向量化运算。空间滞后项必须在中心化之后构造% 中心化之后再乘 W ydm demean_panel(y, X, D, id, time); % 取矩阵下面再拆分 y_c ydm(:, 1); % 第一个输出是 demean 后的 y D_c ydm(:, 3); % 第三个输出是 demean 后的 D X_c ydm(:, 2); % 第二个输出是 demean 后的 X 矩阵 % 生成块对角空间滞后矩阵 Wbig kron(eye(T), W); % 与“时间外层、个体内层”的堆叠顺序一致 Wy Wbig * y_c; % 空间滞后项使用中心化后的 y WD Wbig * D_c; % 政策交叉项使用中心化后的 D WX Wbig * X_c; % 控制变量的空间滞后 % 把变量拼成回归矩阵 Z [D_c, WD, X_c, WX];逻辑说明kron(eye(T), W)生成 NT×NT 的块对角矩阵每个对角块都是同一个 N×N 的 W对应第 1 期到第 T 期。这样做之后Wy的第 i 行代表“第 i 个观测所在时期其空间邻居在该期的加权平均 y”。注意必须先 demean 再乘 W这里 Wy 从构造上就不含固定效应残余。参数说明如果 W 是时变矩阵比如用经济距离权重就不能用kron(eye(T), W)而要做成对每个时期用 W_t 计算Wy(t) W_t * y_c(t)再把各期结果纵向拼起来。经济权重矩阵时变的 SDID 比地理权重难处理得多后面避坑章节会再提。4.2 集中似然 MLE 主循环网格搜索 ρ OLS 闭式解固定效应被中心化之后方程变成 y_c ρ Wy Zδ e。对任意给定的 ρ这个方程是线性回归β 有闭式解因此可以先把 ρ 放在一边对每个候选 ρ 算一次 OLS 残差再代入集中对数似然函数选最大的那个 ρ。这种集中似然法是最稳妥的实现方式中间不需要数值求导也不依赖优化工具箱。function [rho, beta, loglik, se] sdid_mle(y_c, Z, Wy, W, T) % y_c: 中心化后的结果变量 % Z: 中心化后的解释变量矩阵第一列 D第二列 WD % Wy: 中心化后的空间滞后项 % W: 原始 N×N 行标准化权重矩阵 % T: 面板期数 N size(W, 1); n length(y_c); % 预计算 W 的特征值用于快速求 log det(I - rho*W) lam eig(full(W)); rho_grid (-0.98:0.005:0.98); loglik zeros(length(rho_grid), 1); beta_rho zeros(size(Z, 2), length(rho_grid)); resid_rho zeros(length(rho_grid), 1); for i 1:length(rho_grid) r rho_grid(i); % 检查 rho 是否越过奇点 if any(1 - r * lam 0) loglik(i) -Inf; continue; end yr y_c - r * Wy; % 给定 rho 下的 OLS b (Z * Z) \ (Z * yr); e yr - Z * b; sig2 (e * e) / n; % 对数似然- n/2 * (log(2*pi) log(sig2)) T * sum(log(1 - rho*lam)) logdet T * sum(log(1 - r * lam)); loglik(i) -0.5 * n * (log(2 * pi) log(sig2)) logdet; beta_rho(:, i) b; resid_rho(i) sig2; end [~, mi] max(loglik); rho rho_grid(mi); beta beta_rho(:, mi); % 如果想输出标准误需要在最优 rho 下用数值黑塞矩阵或 bootstrap se []; end逻辑说明对数似然里最容易被忽视的是logdet T * sum(log(1 - r * lam))。因为大矩阵 I - ρ Wblock 是 I_T ⊗ (I_N - ρW_N)行列式等于 (I_N - ρW) 行列式的 T 次方所以用 W 的特征值算一次再乘 T 就行。直接对 NT×NT 的稀疏矩阵做log(det(...))N 到 1000、T 到 20 时就已经慢得没法接受。参数说明rho_grid 步长设 0.005 是比较平衡的选择精度要求高可以改成 0.001代价是循环时间拉长约五倍。any(1 - r*lam 0)的判断是防止 ρ 越过奇点比如 W 最大特征值是 1ρ0.99 仍安全但 ρ1 时 log 里出现 0行列式爆炸。如果 W 没有行标准化特征值范围不可控这个检查会直接暴露问题。得到网格最优 ρ 后如果想更高精度可以用fminbnd在最优网格点附近再搜一轮实参是 (r) -nll_m(r)其中 nll_m 用同样的残差逻辑写。4.3 直接效应与间接效应分解β 和 θ 都不是最终答案SDID 的估计系数不能直接当边际效应报告。原因在于 y_c ρWy βD θWD ...变换后 y_c (I - ρW)⁻¹ (βD θWD ...)。某个地区 j 的 D 变化一单位不仅直接改变 j 自己的 y还会通过空间乘数影响邻居 y邻居 y 的变化又反哺回 j。总效应矩阵是S (I - ρW)⁻¹ (β I_N θ W)直接效应 S 对角线的平均值间接效应 S 每行行和减去对角线的平均值总效应等于两者之和。matlab 代码如下。function [direct, indirect, total] sdid_effects(rho, betaD, thetaWD, W) % rho: 空间自回归系数 % betaD: 变量 D 的估计系数 % thetaWD: 变量 WD 的估计系数 % W: N×N 行标准化空间权重矩阵 N size(W, 1); I eye(N); % 空间乘数矩阵 Ainv inv(I - rho * W); % 总效应矩阵 S Ainv * (betaD * I thetaWD * W); direct mean(diag(S)); row_sum sum(S, 2); indirect mean(row_sum - diag(S)); total direct indirect; end逻辑说明S的第 (i,j) 元素表示“j 地区处理状态变化一单位对 i 地区结果产生的累计影响”这个累计已经通过 Ainv 把一轮一轮的间接反馈全部加进来了。diag(S)取的是每个地区对自身的总效应行和减去自身就是对外的溢出。参数说明间接效应的正负号可能和直接效应相反比如一个地区处理挤压了邻居产出direct 为正indirect 为负说明政策有挤出型外溢如果两者同号则属于辐射型外溢。这个分解只对行标准化 W 成立如果你的 W 没有标准化间接效应的大小就没有“邻居加权平均”的解释结果会失真。至于 standard errors我一般用时期块 bootstrap按时间整块重抽样重复跑第 4.1 到 4.2 步得到多组 direct 和 indirect再取标准差。重抽样次数 199 或 399 次够用比直接求解析黑塞矩阵省心很多。5. 亲测避坑五个最容易让 SDID 返工的现场这一章列几个我实际跑 SDID 时踩过或看别人踩过的问题。它们不报红字错误只是让结果在数值上“看起来很对”所以杀伤力比语法错误大得多。每条按现象、原因、解决来写。5.1 空间权重矩阵里的“无邻居”个体行标准化后 NaN 静默传播现象估计结果里 ρ 为 NaN或者 logdet 出现复数往上追查发现 W 里有 NaN。有时候 MATLAB 不报错因为sum(W,2)除出来的 Inf/NaN 会一直传递到特征值和似然函数。原因K 近邻权重矩阵里如果一个城市周边 K 个邻居全是另一个行政区的重复坐标或者经纬度缺失该行距离计算后全为 0行标准化时分母为 0。尤其常见于小岛、飞地、跨境数据。解决标准化之前先执行find(sum(W,2)0)找到零行后按第 3.3 节的方法把自身设成邻居再做标准化。跑完 W 构造代码后用assert(~any(isnan(W(:))))和max(abs(sum(W,2)-1))做双重保险这样后面所有环节默认 W 是干净的。5.2 中心化顺序反过来Wy 用的是未去固定效应的 y现象SDID 跑完ρ 高达 0.9 以上D 的系数符号明显违背直觉而且和普通 DID 差距巨大。你反复检查数据没发现错但总觉得拟合好得不正常。原因写代码时贪图省事先算了Wy Wbig * y再对 Wy 做 demean。因为 demean 投影矩阵 J 与空间权重 W 不可交换加载在未中心化 y 上再做中心化固定效应的时间均值会混进 Wy极大似然会把固定效应误认为空间溢出ρ 被显著高估。解决严格按 4.1 的顺序先把 y、D、X 全部中心化再乘 W。记忆方法就一条先 demean后乘 W。如果非要先乘 W 再 demean需要对方程整体做正交变换比如 Lee-Yu 谱方法那已经超出多数应用场景的需求了不要轻易尝试。这个坑是我的血泪经验检查起来最快的方法是对比Wbig * y_c与demean(Wbig * y)两个向量的相关性如果相关系数明显小于 1说明你的顺序有问题。5.3 行列式在网格边界变成 NaN 或复数ρ 越过了奇点现象loglik 数组在 ρ 接近 ±1 的位置突然变成 NaN或者出现复数网格最优解被顶到边界 0.98 附近似乎模型在暗示“越大越好”。原因W 未行标准化时特征值范围不是 [-1,1]比如如果 W 没有除行和最大特征值可能到 3、5那么 ρ 在 0.5 附近就已经越过奇点。行标准化后特征值最大值是 1但最小值可能是 -1所以 ρ→-1 时也会有奇点。另一个次要原因是 grid 步长太粗恰好落在奇点旁边。解决构造完 W 先跑eig(full(W))观察最大最小特征值。网格上限不要取 0.98 或 0.99 固定值而是取0.99 / max(abs(lam))并确保该值小于 1。同时在似然函数里加上if any(1 - r * lam 0), loglik-Inf; continue; end的判断把越界点直接淘汰。这样即使数据异常也只是 loglik 出现一段 -Inf不会污染最后的选择。5.4 中文注释乱码和数据列名怪字符matlab 环境下最浪费时间的翻车现象拿到的原始代码注释是中文打开后一片乱码或者 CSV 的表头是中文readtable 读进来后列名变成奇怪的制表符导致代码里data.政策没法索引。原因matlab 不同发行版对 UTF-8 和 GBK 的处理策略不统一新版本默认 UTF-8老版本默认本地编码。换电脑、换语言包后代码文件的编码没有跟着转换就会出现乱码。CSV 同理Excel 保存的 CSV 常常是 ANSI 编码。解决拿到代码先做一件事全选、另存为 UTF-8 编码。然后在 matlab 偏好设置里把字符编码改为 UTF-8。如果还乱就直接把注释替换成英文逻辑不受影响。数据列名统一用英文小写加下划线读取时用PreserveVariableNames, true可以避免 matlab 自动把非法字符替换成下划线。这一步不算技术含量但处理的代码包一多这是第一个会让“亲测可用”变成“亲测报错”的地方。5.5 估计结果和普通 DID 一模一样先别高兴很可能是 W 没起作用现象SDID 跑完ρ 的估计值约等于 0 且不显著θ 也不显著β 和普通 DID 完全一致。看起来稳健实际上模型等于没做空间处理。原因最常见的是 W 构造有问题比如mink的索引传反了邻居选成了距离最远的点或者 K 近邻的 K 取 1导致 W 极度稀疏每个地区只连一个邻居空间滞后项几乎没有变异。另一个原因是数据本身无空间自相关但这种情况发生率没那么高。解决先画出莫兰散点图做检验。临时用一个简单脚本计算 y_c 和 Wy 的相关系数如果接近 0说明空间滞后项没有解释力W 选得有问题。再检查 W 每行非零元素个数分布用full(sum(W ~ 0, 2))看一眼正常 K5 时每行应该有 5 个左右非零。如果 W 没问题再看数据里政策虚拟变量的空间分布是不是处理组和对照组在地理上完全混在一起、没有任何空间聚类。最后不是所有题目都适合 SDID有时候普通 DID 的结果就是真的别为了加空间而加空间。6. 不止能跑让 SDID 结果站得住的三件小事代码能出数只是一个开始真正敢写结论前我习惯再做三个验证。第一个是蒙特卡洛自检模拟一个已知 ρ、β、θ 的数据把估计程序跑一遍看能否恢复真实参数。随机生成 y 的过程不难关键还是把第 4 章那套矩阵运算反过来用。% 模拟已知参数的 SDID 数据 rho0 0.5; beta0 1.2; theta0 -0.6; D (rand(N, T) 0.7); D D(:); X randn(N*T, 1); % 按 SDM 公式生成 y注意固定效应和噪声 mu randn(N, 1); mu repmat(mu, T, 1); lambda randn(T, 1); lambda kron(lambda, ones(N, 1)); nu randn(N*T, 1); y (I_NT - rho0 * Wbig) \ (beta0 * D theta0 * (Wbig * D) X mu lambda nu);逻辑说明(I_NT - rho0 * Wbig) \ (...)是模拟内生空间过程的标准写法线性解出来后 y 天然带有空间乘数效应。参数说明蒙特卡洛的样本量不要太小N100、T10 时重复 100 次观察 β 和 θ 的均值是否在真实值附近。如果均值偏出 10% 以上优先检查中心化顺序和 W 标准化。第二个习惯是换 W 做稳健性K4 换到 K8再换反距离阈值权重核心参数的符号和显著性不发生剧烈翻转才能说结果不是一个 K 值选择导致的巧合。第三个习惯是报告效应分解时带上 bootstrap 标准误而不是只报 β 和 θ 的标准误。我自己的流程是拿到一组数据先跑普通 DID 作为基准再跑 SDID两者差距过大就回头查 W 和中心化顺序确认没问题后用蒙特卡洛锁定程序最后写论文时才报告直接效应和间接效应。这套流程走下来虽然慢但基本不会给出一个事后别人复现不了的结论。希望帮到你。本文还有配套的精品资源点击获取