企业官网建设流程全解析

简介：这套代码为空间双重差分模型提供了完整的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 框架。

实际落地时，我一般先跑一个包含 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 排好。最终向量里的顺序是 t=1 期的 N 个观测，接着 t=2 期的 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, K+1, 2); % 多取一个，第 1 个是自身 W = zeros(N); for i = 1:N ne = idx(i, 2:K+1); % 去掉自身 dist_ne = Dk(i, 2:K+1); dist_ne(dist_ne < 1e-6) = 1e-6; % 避免除零 W(i, ne) = 1 ./ dist_ne; % 反距离权重 end % 行标准化放在主函数处理，本函数只返回未标准化矩阵 end

逻辑说明：pdist2计算两两欧氏距离，mink(D, K+1, 2)表示对每一行取最小的 K+1 个值，返回距离矩阵 Dk 和对应列索引 idx。由于每行第一列必然是自身，所以从第 2 列开始截取，就得到 K 个邻居。权重取距离倒数，距离越近影响越大。参数说明：K 的经验范围是 4 到 10，K 太小会让 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 消掉，然后对剩余方程做极大似然估计。这个方案的优点是不需要生成 N+T 个虚拟变量，省内存且收敛快；缺点是对中心化顺序很敏感，顺序错了整个估计结果都是虚的。

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))看一眼，正常 K=5 时每行应该有 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 天然带有空间乘数效应。参数说明：蒙特卡洛的样本量不要太小，N=100、T=10 时重复 100 次，观察 β 和 θ 的均值是否在真实值附近。如果均值偏出 10% 以上，优先检查中心化顺序和 W 标准化。第二个习惯是换 W 做稳健性：K=4 换到 K=8，再换反距离阈值权重，核心参数的符号和显著性不发生剧烈翻转，才能说结果不是一个 K 值选择导致的巧合。第三个习惯是报告效应分解时带上 bootstrap 标准误，而不是只报 β 和 θ 的标准误。

我自己的流程是：拿到一组数据先跑普通 DID 作为基准，再跑 SDID，两者差距过大就回头查 W 和中心化顺序，确认没问题后用蒙特卡洛锁定程序，最后写论文时才报告直接效应和间接效应。这套流程走下来虽然慢，但基本不会给出一个事后别人复现不了的结论。希望帮到你。

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