贸易 + 迁移 + 土地市场：动态空间一体化一般均衡
本页自包含：经济环境与假设 → EK 贸易 + 动态 Bellman 迁移 + 住房供给弹性 → 三市场联合均衡 → 帽子代数动态扩展 → 参数估计 → 福利 → 反事实（贸易自由化 + 迁移联合效应）→ 论文案例。
01 为什么必须联合建模三个市场
在经典贸易模型（Armington、EK、Melitz）中，人口 $L$ 外生、土地市场省略，均衡由贸易平衡与劳动力出清即可封闭。这一简化在跨国贸易问题上近似成立（跨国迁移受签证、语言、文化摩擦极大），但在跨区/跨城市问题上严重失实：
- 贸易成本下降 → 出口行业工资上升 → 工人迁入该地区 → 住房需求上升 → 房价上涨；
- 房价上涨 → 抵消名义工资上升的实际福利 → 抑制进一步迁入；
- 土地所有者（本地居民）从房价上涨中获益，新迁入工人承担"住房成本"；
- 若住房供给弹性低（北京、上海、香港），贸易开放的福利很大一部分被土地所有者"抽走"；若住房供给弹性高（美国 Sun Belt），贸易开放更能转化为实际工资增长。
因此贸易—迁移—土地三者必须同时建模。Caliendo, Parro, Rossi-Hansberg & Sarte (2018, RES；姊妹篇 CP 2015 RES 投入产出版；Caliendo-Dvorkin-Parro 2019 Econometrica 动态版) 的核心贡献就是把这三个市场放进同一个递归系统。
与上一页 RRH 的"通勤分离"不同，本页假设工人住哪就在哪工作，但把跨区迁移内生化。这一选择的现实对应是：在美国，居民在州际间迁移的摩擦主要是"搬家"而非"跨城通勤"；在中国，省际迁移同样是长期居住选择。两个模型的分工是：RRH 回答"城市内部空间结构"（一个都会区内的通勤流），本页回答"跨区域劳动力再配置"（城市之间的人口流）。做中国高铁、都市圈、区域协调发展这类政策评估时，通常需要把两层嵌套：内层 RRH 通勤、外层本页迁移。
贸易成本↓ → 本地市场准入↑ → w↑ → 迁入 L↑ → 住房需求↑ → Q↑ → 实际工资 w/Q ↓ → 迁入放缓。整条链闭环后，福利效应被房价"吸收"一部分。
02 经济环境与假设
模型是动态、多区、多部门的。经济由 $N$ 个地区（如美国 50 州）和 $J$ 个部门组成。时间 $t=0,1,2,\dots$ 离散。核心假设：
- A1：工人可跨区迁移。每期期初可选择留在原地或迁往任意 $i$，迁移成本 $\tau_{ni}\geq 1$（效用损失），idiosyncratic 偏好冲击服从 Fréchet。
- A2：住哪就在哪工作。与 RRH 2017 的"通勤分离"不同，本页工人无通勤，住 $n$ 即工作 $n$。
- A3：住房供给有弹性。住房由土地 + 建筑投入生产，供给弹性 $\varepsilon_n$ 随地区而异。
- A4：商品贸易服从 EK。生产率 Fréchet，贸易份额闭式。
- A5：贴现因子 $\beta\in(0,1)$。迁移决策看当期效用 + 未来延续价值。
03 贸易模块（EK 多部门）
每个地区 $n$、部门 $j$ 的生产率 $z_n^j(\omega)$ 服从 Fréchet：$\Pr(Z_n^j\leq z)=\exp(-T_n^j z^{-\theta_j})$。单位成本：
$\beta_n^j$ 是劳动份额，$\gamma_n^{k,j}$ 是部门 $j$ 对部门 $k$ 的中间投入份额。贸易份额：
04 迁移决策的动态 Bellman
工人在每期期初可选择留在原地区或迁移到任意 $n$。迁移成本 $\tau_{ni}\geq 1$，idiosyncratic 偏好冲击 $\epsilon_{j,n}$ 服从 Fréchet：$G(\epsilon)=\exp(-\epsilon^{-\nu})$。
4.1 Bellman 方程
$U_i(s_t)$ 是 $i$ 地当期间接效用（消费+住房），$\beta$ 是贴现因子。注意：值函数项 $\beta\mathbb{E}_t V_{i,t+1}$ 把预期纳入迁移——修高铁的预期会引发今天的提前迁入。
4.2 迁移概率闭式解 · 完整积分推导（不跳步）
工人从 $n$ 迁到 $i$ 的当期综合收益（含未来延续价值）记为
工人选 $\arg\max_i\{u_{ni}+\epsilon_{ni}\}$，其中 idiosyncratic 偏好冲击 $\epsilon_{ni}$ 独立同 Type-I 极值（Gumbel）分布：
Step 1 · 事件概率。工人选 $i$ 当且仅当 $u_{ni}+\epsilon_{ni}\ge u_{nk}+\epsilon_{nk}\;\forall k$。以 $\epsilon_{ni}=v$ 为条件（密度 $g(v)$），要求其余 $\epsilon_{nk}\le v-(u_{ni}-u_{nk})$：
Step 2 · 代入 CDF。把 $G(v-u_{ni}+u_{nk})=\exp(-e^{-(v-u_{ni}+u_{nk})/\nu})=\exp(-e^{-v/\nu}e^{u_{ni}/\nu}e^{-u_{nk}/\nu})$ 连乘，并把 $k=i$ 自己那一项也并入指数：
Step 3 · 变量替换 $z=e^{-v/\nu}$。由 $z=e^{-v/\nu}$ 得 $v=-\nu\ln z$，微分 $dz=-(1/\nu)e^{-v/\nu}dv$，即 $\frac{1}{\nu}e^{-v/\nu}dv=-dz$。换上下限：$v:-\infty\to+\infty$ 对应 $z:+\infty\to 0$：
其中用了 $\int_0^\infty e^{-az}dz=1/a$（$a>0$）。
4.3 人口运动定律与稳态
给定初始分布 $\{L_{n,0}\}$，这条方程完整描述人口动态演化。稳态：$L_{n,t+1}=L_{n,t}=L_n^*$，即 $L^*$ 是迁移矩阵 $M=[m_{ni}]$ 的主特征向量：
05 土地市场与住房供给弹性
5.1 住房生产
5.2 住房供给弹性 · 正式 reduced-form 方程
住房供给方最大化利润 $Q_n H_n - r_n L_n^{\text{land}} - w_n I_n^H$。一阶条件给出 $I_n^H$ 与 $L_n^{\text{land}}$ 的需求，结合外生土地供给，可推出房价对人口的 reduced-form 弹性。对数线性化后：
5.3 住房市场出清
注意：与 RRH 2017 不同，这里工人住在哪就在哪工作（无通勤），所以 $n$ 地工资 $w_n$ 就是 $n$ 地居住工人的工资。
06 三市场联合均衡
一个动态空间一般均衡由 $\{w_{n,t}, P_{n,t}^j, Q_{n,t}, L_{n,t}, \pi_{ni,t}^j, m_{ni,t}, V_{n,t}\}$ 组成，满足：
- 企业最大化（$c_n^j, \pi_{ni}^j$）；
- 工人迁移优化（Bellman + $m_{ni,t}$）；
- 产品市场出清（部门 $j$ 贸易平衡）；
- 劳动力市场出清：$w_n L_n = \sum_j \beta_n^j [\sum_i \pi_{in}^j X_i^j]$；
- 土地市场出清：$Q_n H_n=(1-\beta)w_n L_n$；
- 人口运动定律：$L_{i,t+1}=\sum_n m_{ni,t}L_{n,t}$。
求解采用动态迁移 + 静态贸易均衡的联合迭代：外循环同时更新 $\{w_n, Q_n, L_n\}$，内循环在给定 $w$ 下求住房供给曲线与迁移概率，再用不动点迭代稳态人口。
数值实现上，静态贸易均衡对给定工资是闭式的（贸易份额、价格指数、住房出清都可以直接算），难点在动态迁移的稳态人口分布。标准做法是先向后迭代值函数（从 Terminal 期向前收敛 Bellman），再向前迭代人口运动定律直到时间不变。稳态反事实则可以用帽子代数跳过转移动态，直接比较两个稳态；做动态反事实（如"关税在 $T=5$ 突然下降"）则必须完整模拟人口从初始分布到新稳态的过渡路径，通常需要 50–100 期迭代。
07 帽子代数的动态扩展
反事实用"帽子代数"避免水平值校准，只需基准均衡份额。完整七方程如下：
(i) 单位成本帽子（投入产出链）：
$$\hat c_n^j=(\hat w_n)^{\beta_n^j}\prod_k(\hat P_n^k)^{\gamma_n^{k,j}}$$(ii) 贸易份额帽子：
$$\hat\pi_{ni}^j = \frac{\pi_{ni}^j (\hat c_i^j \hat d_{ni}^j)^{-\theta_j}}{\sum_k \pi_{nk}^j (\hat c_k^j \hat d_{nk}^j)^{-\theta_j}}$$(iii) 价格指数帽子：
$$\hat P_n^j = \left[\sum_i \pi_{ni}^j (\hat c_i^j \hat d_{ni}^j)^{-\theta_j}\right]^{-1/\theta_j}$$(iv) 住房出清帽子（含供给弹性 $\varepsilon_n$）：
$$\hat Q_n = \hat L_n^{1/\varepsilon_n}$$(v) 间接效用帽子：
$$\hat U_n=\frac{\hat w_n}{(\hat P_n)^\beta(\hat Q_n)^{1-\beta}}$$(vi) 迁移概率帽子（含未来值函数项）：
$$\hat m_{ni}=\frac{\exp\!\left[\frac{1}{\nu}(\log\hat U_i-\log\hat\tau_{ni}+\beta\hat V_i)\right]}{\sum_k m_{nk}\exp\!\left[\frac{1}{\nu}(\log\hat U_k-\log\hat\tau_{nk}+\beta\hat V_k)\right]}$$(vii) 人口运动帽子（动态）：
$$\hat L_{i,t+1}=\frac{\sum_n m_{ni,t}L_{n,t}\hat m_{ni,t}}{\sum_n m_{ni,t}L_{n,t}}$$稳态反事实可用 (i)–(vii) 联立迭代 $\hat w,\hat L,\hat Q$ 直到收敛；动态反事实（如 $T=5$ 关税突然下降）则要在 (vi)(vii) 之间向前模拟 50–100 期过渡路径，每期都重解静态均衡。
当迁移成本 $\tau_{ni}$ 本身变化时需知道基准 $m_{ni}$ 的水平值；当模型有状态变量（累积资本）时，帽子代数要扩展为动态递归，不能一次性静态求解；当 $\varepsilon_n$ 地区异质时，(iv) 必须逐地区代入，不能用全国平均弹性。
08 参数估计
8.1 贸易弹性 $\theta_j$
按部门估计：用部门 $j$ 双边贸易流 + 关税变化做 PPML 回归（见 03 页），系数即 $-\theta_j$。典型值：制造业 $\theta_j\approx 4\sim 5$，服务业 $\theta_j$ 较大（~10）。
8.2 迁移弹性 $\nu$
对迁移概率取对数：$\ln(m_{ni}/m_{nn}) = \frac{1}{\nu}(\log U_i - \log\tau_{ni} + \beta V_i) - \frac{1}{\nu}\log\tau_{nn}$。用人口普查/年度迁移流数据 + 观测实际工资与房价，加 origin 固定效应吸收 $\tau_{nn}$，即可 OLS/IV 估计 $1/\nu$。典型 $\nu\approx 1.5\sim 3$（Kennan-Walker 2011; Monar-Redding 2022）。
8.3 住房供给弹性 $\varepsilon_n$
用 Saiz (2010) 方法：城市级面板回归 $\Delta\log Q_{n,t}$ 对 $\Delta\log L_{n,t}$，工具变量用地理约束（不适合建楼的土地比例）与建筑管制强度（Wharton Residential Land Use Regulation Index）。典型 $\varepsilon_n$：曼哈顿 ~0.5，休斯顿 ~3。
8.4 数据要求与具体来源
| 数据 | 国际来源 | 中国来源 | 用途 |
|---|---|---|---|
| 双边贸易流（分部门） | CEPII BACI（HS6 分档）、UN Comtrade、IMF DOTS、OECD-ICIO / WIOD 投入产出表 | 中国海关数据库（2000–2016 企业级进出口）；中国 31 省投入产出表（国家统计局）；各省间投入产出表 | 估计 $\theta_j$、构造 $\pi_{ni}^j$ |
| GDP/工资/就业 | PWT 9.1+（Penn World Table）、World Bank WDI、ILO | 《中国统计年鉴》、《中国城市统计年鉴》、CEIC、EPS 数据平台 | $w_n$、$L_n$、部门就业 |
| 双边迁移流 $m_{ni}$ | OECD International Migration Database、US Census ACS 1-year migration flow、IRS SOI 税号迁移 | 人口普查长表（2000/2010/2020 户籍地-常住地匹配）；流动人口动态监测调查 CMDS（卫健委，个体级面板） | 估计 $\nu$、构造 $m_{ni}$ |
| 房价/租金 $Q_n$ | BIS Residential Property Price、OECD House Prices、Zillow、CPPI | 国家统计局 70 城新建/二手住宅价格指数；中国房价行情网 CREIS；Wind/CEIC 城市房价面板 | 住房出清、$\varepsilon_n$ 估计 |
| 住房供给弹性 $\varepsilon_n$ | Saiz (2010) 城市级估计；Wharton WRLURI 管制指数；USGS 地形坡度 | 城市地形坡度（CGIAR SRTM 卫星）；土地利用现状变更调查；城市建成区面积（《城市建设统计年鉴》） | reduced form $\hat Q_n=\hat L_n^{1/\varepsilon_n}$ |
| 双边贸易成本 $d_{ni}^j$ | CEPII Gravity（距离/接壤/共同语言）；关税 WITS/TAO | 铁路/公路最短里程（OpenStreetMap+OSRM）；省间铁路货运量；关税税则 | EK 贸易份额 |
| 部门间投入产出 | WIOD 2016、OECD-ICIO、GTAP | 国家统计局全国投入产出表（每 5 年）；省间 IO 表（各省统计局） | $c_n^j$ 中 $\gamma_n^{k,j}$ |
实操提醒：(i) 中国省际贸易流没有官方连续统计，常用"引力模型反推"或"省间铁路货运量+公路客货周转"近似；(ii) CMDS 是唯一能同时识别迁移成本与自选择的个体面板，但 2018 年后公开版本缩减；(iii) $\varepsilon_n$ 中国版没有现成 Saiz 估计，需用卫星坡度 + 土地利用管制指数自建；(iv) 动态求解要求每期都有一套均衡，因此数据频率最好是年频，避免用季度数据引入季节噪声。
09 福利分析
9.1 工人福利
间接效用 $U_n = w_n/(P_n^\beta Q_n^{1-\beta})$。在空间均衡中，跨区效用（含迁移选项价值）由值函数 $V_n$ 刻画。平均福利 $\bar U = \sum_n L_n U_n / \bar L$。
9.2 迁移的福利效应
迁移自由化的福利效应来自三个渠道：
- (a) 选址效率收益：工人从低生产率地迁到高生产率地，总产出↑；
- (b) 房价反馈：流入地房价↑、流出地↓，部分抵消名义工资差异；
- (c) 跨期平滑：动态 Bellman 中工人可"提前搬家"以预期未来收益，福利反应比静态模型更平滑。
9.3 关键发现（CPRS 2018）
忽略土地市场的模型会系统性高估贸易开放的工人福利，因为没考虑房价上涨对实际工资的侵蚀。在美国州级数据中，土地市场反馈吸收了约 20–30% 的名义贸易收益；住房供给弹性低的州（如加州沿海）吸收比例更高。
9.4 分配维度：工人 vs 土地所有者
总福利 $\bar U$ 的平均化会掩盖巨大的分配差异。贸易自由化后，本地有房的原住民（土地所有者）从房价上涨中获益，新迁入工人则承担高房价。做政策评估时，必须同时报告"总福利变化"与"按居住地/按存量房主/按新移民"的分组福利变化，否则会得出"政策对所有人都好"的错误结论。CPRS (2018) 发现，美国制造业冲击后，铁锈地带（Rust Belt）的房价下跌使当地存量房主财富缩水，而工人可通过迁移到 Sun Belt 对冲——但住房供给弹性低的地区工人无法充分调整，承担了大部分调整成本。
10 Python 求解器与反事实
本节的反事实把帽代数扩展到 $\hat w_i$ + $\hat L_i$ + 地租/房价三个市场。通用帽代数骨架、阻尼 / Anderson 加速、从 EK → RRH → 三市场的嵌套迁移差异 → 先学 基础知识库 · 帽代数反事实求解
下面是静态简化版"贸易+迁移+土地"联合求解器（单部门、无中间投入，保留核心反馈），并附反事实代码。
# ============================================================
# 贸易 + 迁移 + 土地 一体化一般均衡求解（静态简化版）
# 框架参考: Caliendo-Parro-Rossi-Hansberg-Sarte (2018)
# ============================================================
import numpy as np
def solve_tml(N=6, theta=4.0, nu=2.0, beta=0.75,
eland=2.0, damp=0.5, tol=1e-7,
maxiter=3000, seed=0, d_shock=None, tau_shock=None):
"""
静态简化版：EK 贸易 + Frechet 迁移 + 住房供给弹性
d_shock, tau_shock: 可选 (N,N) 反事实冲击
"""
rng = np.random.default_rng(seed)
T = rng.uniform(0.8, 1.2, N)
H0 = rng.uniform(0.8, 1.2, N)
Lbar = 1.0
dist = rng.uniform(0.5, 2.0, (N, N))
d = np.exp(0.3 * dist)
np.fill_diagonal(d, 1.0)
tau = np.exp(0.5 * dist)
np.fill_diagonal(tau, 1.0)
if d_shock is not None:
d = d * d_shock
if tau_shock is not None:
tau = tau * tau_shock
w = np.ones(N); L = np.ones(N)/N; Q = np.ones(N)
for it in range(maxiter):
# 1. 价格指数 + 贸易份额
A = T[None, :] * (w[None, :] * d) ** (-theta)
P = np.sum(A, axis=1) ** (-1.0/theta)
pi = A / np.sum(A, axis=1, keepdims=True)
# 2. 住房出清：Q^(1+1/eland) = (1-beta) w L / H0
Q_new = ((1-beta) * w * L / H0) ** (1.0/(1.0+1.0/eland))
# 3. 间接效用
U = w / (P**beta * Q_new**(1-beta))
# 4. 迁移概率 m_ni
V = np.exp((1.0/nu)*(np.log(U)[None, :] - np.log(tau)))
m = V / np.sum(V, axis=1, keepdims=True)
# 5. 稳态人口不动点
L_new = L.copy()
for _ in range(50):
L_next = np.sum(m * L_new[:, None], axis=0)
L_next = L_next / L_next.sum() * Lbar
if np.max(np.abs(L_next - L_new)) < 1e-10:
L_new = L_next; break
L_new = 0.5*L_next + 0.5*L_new
# 6. 工资方程（贸易平衡）
sales_i = np.sum(pi * (beta*w*L_new)[:, None], axis=0)
w_new = sales_i / L_new
w_new = w_new / w_new.mean()
Q = damp*Q_new + (1-damp)*Q
w = damp*w_new + (1-damp)*w
L = L_new
err = np.max(np.abs(w - w_new)) + np.max(np.abs(Q - Q_new))
if err < tol:
print(f"收敛于第 {it} 轮, err={err:.2e}"); break
return dict(w=w, Q=Q, L=L, P=P, pi=pi, m=m, U=U)
if __name__ == "__main__":
# 基准
base = solve_tml(N=5)
# 反事实：所有贸易成本与迁移成本 -20%
shock = 0.8 * np.ones((5, 5))
np.fill_diagonal(shock, 1.0)
cf = solve_tml(N=5, d_shock=shock, tau_shock=shock)
print("基准福利 U:", np.round(base["U"], 4))
print("反事实福利 U':", np.round(cf["U"], 4))
print("福利变化 U'/U:", np.round(cf["U"]/base["U"], 4))
# 直觉：贸易+迁移自由化 → 平均福利↑，但流入地房价↑
10.2 完整动态求解器：值函数向后迭代 + 人口向前运动 + 稳态主特征向量
静态版只算稳态，不能回答"政策冲击后多少年过渡完、过渡期福利路径如何"。下面是 Caliendo-Dvorkin-Parro (2019) 式的完整动态求解器：
# ============================================================
# 动态"贸易+迁移+土地"联合求解器
# 框架参考: Caliendo-Dvorkin-Parro (2019, Econometrica)
# 流程:
# (A) 静态内核: 给定 L_t, 解当期均衡 (w_t, P_t, Q_t, U_t)
# (B) 向后迭代: 从 terminal T 倒推值函数 V_t
# (C) 向前迭代: 由 m_{ni,t} 推人口 L_{t+1}
# (D) 稳态: L* = M^T L* (主特征向量 / 反复迭代)
# ============================================================
import numpy as np
# ---------- (A) 静态内核: 给定 L, 解当期 w,P,Q,U ----------
def static_equilibrium(L, T, d, theta, beta, eland, H0,
w_init=None, tol=1e-8, maxiter=2000, damp=0.4):
"""给定人口分布 L, 不动点迭代 w 与 Q。返回当期 U。"""
N = len(L)
w = np.ones(N) if w_init is None else w_init.copy()
Q = np.ones(N)
for it in range(maxiter):
A = T[None, :] * (w[None, :] * d) ** (-theta)
P = np.sum(A, axis=1) ** (-1.0 / theta)
pi = A / np.sum(A, axis=1, keepdims=True)
# 住房供给弹性 reduced form: Q^(1+1/eland) = (1-beta) w L / H0
Q_new = ((1 - beta) * w * L / H0) ** (1.0 / (1.0 + 1.0 / eland))
U = w / (P ** beta * Q_new ** (1 - beta))
# 工资方程(贸易平衡): w_i L_i = sum_n pi_ni * (beta * w_n L_n)
sales_i = np.sum(pi * (beta * w * L)[:, None], axis=0)
w_new = sales_i / L
w_new = w_new / w_new.mean()
err = np.max(np.abs(w_new - w)) + np.max(np.abs(Q_new - Q))
w = damp * w_new + (1 - damp) * w
Q = damp * Q_new + (1 - damp) * Q
if err < tol:
break
return dict(w=w, P=P, Q=Q_new, U=U, pi=pi)
# ---------- (B) 向后迭代值函数 V_t ----------
def backward_iteration(U_path, tau, beta, nu, T_horiz):
"""
U_path[t] = N-向量, 各期间接效用(由静态内核在候选政策路径下解出)
V_T = U_T (terminal value)
V_n,t = nu * log sum_i exp((log U_i,t - log tau_ni + beta V_i,t+1)/nu)
(Type-I 极值的 log-sum-exp 闭式期望)
"""
N = len(U_path[0])
V = [None] * (T_horiz + 1)
V[T_horiz] = np.log(U_path[T_horiz]) # terminal
for t in range(T_horiz - 1, -1, -1):
# u_ni = log U_i,t - log tau_ni + beta V_i,t+1 (N x N)
u = np.log(U_path[t])[None, :] - np.log(tau) + beta * V[t + 1][None, :]
# m_ni = exp(u_ni/nu) / sum_k exp(u_nk/nu), 行归一化
u_scaled = u / nu
u_max = u_scaled.max(axis=1, keepdims=True)
logsum = u_max[:, 0] + np.log(np.sum(np.exp(u_scaled - u_max), axis=1))
V[t] = nu * logsum # = E[max_i u_ni + eps_ni]
return V
# ---------- (C) 向前人口运动 ----------
def forward_population(V, U_path, tau, nu, L0, T_horiz):
L_path = [None] * (T_horiz + 1)
L_path[0] = L0.copy()
for t in range(T_horiz):
u = (np.log(U_path[t])[None, :] - np.log(tau)
+ V[t + 1][None, :]) / nu
u -= u.max(axis=1, keepdims=True)
m = np.exp(u)
m = m / m.sum(axis=1, keepdims=True) # m_ni: 行 n 列 i
L_next = np.sum(m * L_path[t][:, None], axis=0)
L_path[t + 1] = L_next / L_next.sum() * L0.sum()
return L_path, m
# ---------- (D) 稳态主特征向量 ----------
def steady_state_population(m, Lbar=1.0, tol=1e-12, maxiter=100000):
"""反复迭代 L <- m^T L 直到收敛 (Perron-Frobenius)。"""
L = np.ones(len(m)) / len(m)
for _ in range(maxiter):
L_next = m.T @ L
L_next = L_next / L_next.sum() * Lbar
if np.max(np.abs(L_next - L)) < tol:
return L_next
L = L_next
return L
# ---------- 主程序: 一次反事实 = 静态均衡路径 + 动态过渡 ----------
def solve_dynamic(N=5, T_horiz=60, theta=4.0, nu=2.0, beta_disc=0.95,
beta=0.75, eland=2.0, seed=0, tau_shock=None):
rng = np.random.default_rng(seed)
T_prod = rng.uniform(0.8, 1.2, N)
H0 = rng.uniform(0.8, 1.2, N)
dist = rng.uniform(0.5, 2.0, (N, N))
d = np.exp(0.3 * dist); np.fill_diagonal(d, 1.0)
tau = np.exp(0.5 * dist); np.fill_diagonal(tau, 1.0)
if tau_shock is not None:
tau = tau * tau_shock
L0 = np.ones(N) / N # 初始人口
# 1) 先在初始 L 下算静态均衡, 得 U_path (简化: 假设政策路径不变)
U_path = [None] * (T_horiz + 1)
L_path = [None] * (T_horiz + 1)
L_path[0] = L0.copy()
w_guess = np.ones(N)
for t in range(T_horiz + 1):
eq = static_equilibrium(L_path[t], T_prod, d, theta,
beta, eland, H0, w_init=w_guess)
U_path[t] = eq["U"]; w_guess = eq["w"]
if t < T_horiz:
# 先用静态 m 推一步人口, 让 L_path 与 U_path 一致
u = np.log(U_path[t])[None, :] - np.log(tau)
m_s = np.exp(u / nu); m_s /= m_s.sum(axis=1, keepdims=True)
L_next = np.sum(m_s * L_path[t][:, None], axis=0)
L_path[t + 1] = L_next / L_next.sum()
# 2) 向后迭代 V
V = backward_iteration(U_path, tau, beta_disc, nu, T_horiz)
# 3) 用含 V 的 m 重新向前推一次(含预期项)
L_path, m_final = forward_population(V, U_path, tau, nu, L0, T_horiz)
# 4) 稳态: m_final 的主特征向量
L_ss = steady_state_population(m_final)
return dict(U_path=U_path, L_path=L_path, V=V,
m=m_final, L_ss=L_ss, L_path_final=L_path)
if __name__ == "__main__":
sol = solve_dynamic(N=4, T_horiz=40)
print("t=0 人口分布:", np.round(sol["L_path"][0], 3))
print("t=20 人口分布:", np.round(sol["L_path"][20], 3))
print("稳态人口(主特征向量):", np.round(sol["L_ss"], 3))
print("收敛检验 L_ss vs L_path[-1]:",
np.allclose(sol["L_ss"], sol["L_path"][-1], atol=1e-3))
值函数 $V_{i,t+1}$ 出现在 $t$ 期迁移概率里，而 $V_{i,t+1}$ 又依赖 $t+1$ 期及以后的 $U$。这是典型的动态规划结构：从 terminal 期 $V_T=\log U_T$ 向前倒推（backward induction），得到整条 $\{V_t\}$ 路径后，再用 $t=0$ 初始人口向前推 $\{L_{t+1}\}$。这就是 Caliendo-Dvorkin-Parro (2019) 的标准算法，也是 DSGE 中"Hall/Judd"方案在空间模型上的对应。
10.1 反事实典型结果
| 市场机制 | 对 $\hat U$ 贡献 | 直觉 |
|---|---|---|
| 贸易成本↓ | $+$ | 更多品种、更低 $P_n$ |
| 迁移成本↓ | $+$ | 工人可去更好的地方 |
| 人口重新分布 | $+/-$ | 高工资地涌入 → 房价↑ → 部分抵消 |
| 土地供给弹性 $\varepsilon_n$ | 弹性越大 $+$ 越多 | Sun Belt 获益多于 Manhattan |
政策含义上，把"贸易自由化"与"迁移自由化"放在一起反事实，可以回答一个经典问题：如果中国劳动力可以自由跨区流动，制造业出口扩张的福利会如何分布？CPRS 类模型给出的答案是：贸易自由化的名义收益大部分被沿海城市的房价上涨吸收，工人福利增长远小于名义工资增长；如果同时放松落户限制（降低 $\tau_{ni}$），迁移成本下降会让更多工人向沿海集聚，虽然短期加剧沿海房价压力，但长期通过竞争效应压低租金，并通过区域间工资均等化提升整体效率。区域政策评估的关键，是把"贸易成本"与"迁移成本"作为两个独立政策杠杆，分别模拟。
11 论文案例与常见错误
11.1 论文案例
难度与前置
| 维度 | 要求 |
|---|---|
| 数学 | 动态规划/Bellman 方程、Type-I 极值 logit 闭式、Perron-Frobenius 主特征向量、对数线性化 |
| 计量 | 多部门引力 PPML、迁移流离散选择 MLE、IV（地理工具变量） |
| 编程 | numpy 向量化、向后迭代+向前迭代、log-sum-exp 数值稳定技巧；建议先跑通本目录 06 页 RRH 静态求解器 |
| 前置阅读 | 本目录 02 EK → 04 福利分解 → 05 帽代数 → 06 RRH 静态 → 本页动态扩展 |
| 工作量估计 | 单部门静态版 ≈ 1 周；加 Bellman 动态路径 ≈ 2–3 周；多部门+IO 投入产出 ≈ 1–2 个月 |
| 常见坑 | 见下方 5 条；最坑的是用静态迁移概率代替动态 Bellman 项 |
11.2 常见错误（5 条）
动态模型中迁移概率含未来值函数项 $\beta\mathbb{E}V$。用静态 Frechet 概率代替，会高估短期迁移反应、低估长期调整。
直接设 $Q_n=1$ 或 $H_n$ 固定，会把所有贸易收益转化为名义工资，完全忽略房价反馈，结论严重偏误。
从北京到鹤岗与从鹤岗到北京的迁移成本完全不同。对称化 $\tau_{ni}=\tau_{in}$ 会把"退出成本"与"进入成本"混淆，估计出错误稳态分布。
单部门模型中制造业只用劳动。现实中部门 $j$ 成本依赖部门 $k$ 价格 $P_n^k$。忽略 IO 会低估跨部门溢出（如能源价格↑→所有部门成本↑）。
两者都是 Freéchet 形状参数但识别来源完全不同（迁移流 vs 贸易流），数值也不同（$\nu\approx 2$, $\theta\approx 5$）。混用会导致反事实福利偏估 2–3 倍。
11.3 进阶资料
- Ossa, R. (2016), "Quantitative Models of Commercial Policy", Handbook of Commercial Policy。
- 范剑勇等（2019）《量化空间经济学在中国的应用》综述。
- Reddding (2020), "Workhorse Tools in QSE"。
∑ 方程总清单 · Equation Summary
本模型共 12 个方程，按推导顺序编号如下。
| 编号 | 名称 | 公式 | 所在节 |
|---|---|---|---|
| (TM1) | 多部门单位成本 | $c_n^j=\Upsilon_n^j(w_n)^{\alpha_n^j}\prod_k(P_n^k)^{\gamma_n^{kj}}$ | 03 |
| (TM2) | ★ 多部门贸易份额 | $\pi_{ni}^j=T_i^j(c_i^j\tau_{ni}^j)^{-\theta_j}/\sum_k T_k^j(c_k^j\tau_{nk}^j)^{-\theta_j}$ | 03 |
| (TM3) | 多部门价格指数 | $P_n^j=\gamma_j[\sum_k T_k^j(c_k^j\tau_{nk}^j)^{-\theta_j}]^{-1/\theta_j}$ | 03 |
| (TM4) | ★ Bellman 方程 | $V_{n,t}=\max_i\{\log U_{i,t}-\log\tau_{ni}+\epsilon_{ni,t}+\beta\mathbb E_t V_{i,t+1}\}$ | 04 |
| (TM5) | ★ 迁移概率（闭式） | $m_{ni}=\exp[u_{ni}/\nu]/\sum_k\exp[u_{nk}/\nu]$ | 04 |
| (TM6) | 人口运动定律 | $L_{i,t+1}=\sum_n m_{ni,t}L_{n,t}$ | 04 |
| (TM7) | 稳态人口 | $L^*=M^\top L^*$（主特征向量） | 04 |
| (TM8) | 住房生产 | $H_n=A_n^H(L_n^{land})^{\phi_n}(I_n^H)^{1-\phi_n}$ | 05 |
| (TM9) | ★ 住房供给弹性 | $\hat Q_n=\hat L_n^{1/\varepsilon_n}$ | 05 |
| (TM10) | 住房出清 | $Q_n H_n=(1-\beta)w_n L_n$ | 05 |
| (TM11) | ★ 三市场联合均衡 | 贸易出清 + 迁移稳态 + 住房出清（联立不动点） | 06 |
| (TM12) | 帽代数方程组 | $\hat P_n^j,\hat\pi_{ni}^j,\hat c_i^j,\hat w_n,\hat L_n,\hat Q_n$ | 07 |