Caliendo-Parro (2015)：多国多部门量化贸易模型与投入产出关联
把单部门 EK 升级为 $N$ 国 $J$ 部门、部门间通过投入产出（IO）矩阵相互需求、部门特定贸易弹性 $\theta_j$、从价关税 $\tau$ 与关税收入返还的可量化一般均衡。本页从 Fréchet 分布逐部门复制，到嵌套工资不动点求解，再到 ACR 福利分解与帽子代数反事实，全程自包含。
numpy 的广播与高维张量求和。numpy。本页求解器是"外层工资不动点、内层价格/支出不动点"的双层迭代，比单部门 EK 多一层 IO 耦合。01 经济环境与核心假设
Caliendo & Parro (2015, AER) 把 Eaton-Kortum (2002) 从"单一部门"推广到多国 $N$ × 多部门 $J$，并允许部门之间通过投入产出（IO）矩阵相互需求。这是今天所有量化贸易政策评估（NAFTA、RCEP、中美贸易摩擦、关税）的标准工作母机。
经济环境的四条骨架假设（每条附设定理由 / 参数含义 / 经济直觉）：
- 假设 A1：$N$ 国、$J$ 部门；劳动在一国内部门间可流动、跨国不可流动。 理由：部门间流动保证一国只有一个工资 $w_i$（不论在哪个部门干活），跨国不可流动则工资在国家间内生差异。参数含义：$L_i$ 是 $i$ 国总劳动力；$w_i$ 是 $i$ 国统一工资。直觉：这正是美国工人可以跨行业跳槽、却很难搬到墨西哥的现实。
- 假设 A2：每个部门 $j$ 的生产率服从部门特定的 Fréchet 分布 $F_i^j(z)=e^{-T_i^j z^{-\theta_j}}$。 理由：与 02-EK 相同的极值闭包性质给出闭式份额；但形状参数现在是部门特定 $\theta_j$。参数含义：$T_i^j$ 是 $i$ 国 $j$ 部门绝对生产率；$\theta_j$ 是 $j$ 部门的比较优势离散度，同时就是 $j$ 部门的贸易弹性。直觉：化工品技术差异大（$\theta$ 小、弹性高），农产品差异小（$\theta$ 大、弹性低），不能用同一个 $\theta$。
- 假设 A3：部门间投入产出关联。部门 $j$ 的生产不仅用劳动，还用所有部门 $k$ 的中间品，份额为 $a_{kj}$（中间品篮子内归一，$\sum_k a_{kj}=1$）。 理由：现实中钢材厂买电力、汽车厂买钢材。参数含义：$\beta_j$ 是劳动（增加值）份额，$1-\beta_j$ 是中间品份额；$a_{kj}$ 来自 IO 表。直觉：关税对钢材提价，会沿 IO 链条传到汽车成本——这是单部门 EK 完全漏掉的渠道。
- 假设 A4：冰山成本 $\tau_{ni}^j\ge 1$ 与从价关税 $t_{ni}^j\ge 0$，关税收入一次性返还给家庭。 理由：$\tau$ 概括运费/距离/非关税壁垒；关税是政策变量，$t_{nn}^j=0$。参数含义：到岸价 = $c_i^j \tau_{ni}^j (1+t_{ni}^j)$。直觉：关税是"把钱从消费者口袋转到政府再发回家庭"的转移，不直接破坏效率但改变相对价格。
单部门 EK 里单位成本就是 $w_i$。多部门 CP 里单位成本 $c_i^j$ 是工资与所有部门价格指数的 Cobb-Douglas 组合。这意味着：任何一个部门 $k$ 的价格下降，都会通过中间品渠道降低所有使用 $k$ 中间品的部门 $j$ 的成本——价格在部门网络里全局耦合。这就是为什么必须用"内层价格不动点 + 外层工资不动点"的双层算法。
02 部门单位成本：工资 + 中间品
部门 $j$ 的技术是 Cobb-Douglas：1 单位产出需要 $\beta_j$ 单位"劳动"与 $(1-\beta_j)$ 单位"中间品复合包"，中间品复合包又是各部门价格指数 $P_i^k$ 的 Cobb-Douglas，权重 $a_{kj}$。单位成本（价格等于边际成本，CRS+完全竞争）：
设定理由：Cobb-Douglas 使"份额 = 指数"，单位成本取 log 后变成线性组合，数值稳定。参数含义：$\beta_j$ 越大，部门 $j$ 越"劳动密集"，成本越靠工资；$1-\beta_j$ 是中间品密集度，制造业通常 $\beta_j\approx 0.3$，服务业 $\beta_j\approx 0.6$。经济直觉：$a_{kj}$ 列就是 IO 表中"部门 $j$ 每 1 元中间品支出里 $k$ 部门占多少"。
两边取 log（数值上必须这么做，避免大幂次下溢）：
当 $\tau$、$\theta$ 较大时，$(c_i^j\tau_{ni}^j(1+t_{ni}^j))^{-\theta_j}$ 可能 $<10^{-300}$ 直接下溢成 0。所有 CP 实现都在 log 空间做幂次，最后再 np.exp。本页代码严格这么写。
03 部门价格指数与贸易份额
在部门 $j$ 内，$n$ 国消费者面对所有 $i$ 国的到岸价 $c_i^j\tau_{ni}^j(1+t_{ni}^j)$。与 02-EK 完全相同的 Fréchet 闭包论证（此处不重复积分，只给结果），得到部门 $j$ 的"可竞争成本"总和：
推导不跳步说明：把 02-EK 第 (EK6)–(EK7) 的"最低价分布"逐部门重做一遍——$n$ 国从 $i$ 国买到部门 $j$ 品种的概率，正比于 $T_i^j(c_i^j\tau_{ni}^j(1+t_{ni}^j))^{-\theta_j}$，对 $i$ 加总即分母。价格指数只差一个常数 $\gamma$（含 $\sigma$），反事实里被帽子消掉，故略去。
部门贸易份额：
参数含义：$\pi_{ni}^j$ 是 $n$ 国在部门 $j$ 上花的每一元钱里，从 $i$ 国买的比例（含关税的到岸支出份额）。经济直觉：分子是 $i$ 国在 $j$ 部门的"竞争力得分"（生产率高、工资低、成本低、关税低则大），分母是所有来源国的得分加总。注意 $\sum_i\pi_{ni}^j=1$（行和为 1）。
CP 把关税写成 (1+t) 乘进到岸成本里，与冰山 $\tau$ 并列。差别在于：关税是转移支付（钱留给 $n$ 国政府），冰山是真实损耗。所以在市场出清时，$i$ 国真正收到的是 $CIF$ 价值 $\pi_{ni}^j X_n^j/(1+t_{ni}^j)$，而不是含关税的 $\pi_{ni}^j X_n^j$。这一步最容易错（见常见错误 3）。
04 支出、关税收入与市场出清
部门 $j$ 的总支出 $X_n^j$ 有两个去向：(1) 家庭最终消费；(2) 各部门生产用的中间品需求。设家庭把可支配收入按固定份额 $b_j$ 花在部门 $j$（$\sum_j b_j=1$），则：
逐项解读：右边第一项是最终消费，$w_n L_n+TR_n$ 是家庭可支配收入（劳动收入 + 关税返还）；第二项是中间品需求——部门 $k$ 花在中间品上的总额是 $(1-\beta_k)X_n^k$，其中比例 $a_{jk}$ 用于采购部门 $j$ 的产品，对 $k$ 加总。
关税收入 $TR_n$（政府一次性返还）：
工资市场出清：$i$ 国劳动收入 = 所有国家 $n$ 在所有部门 $j$ 从 $i$ 国购买的 CIF 价值之和：
Walras 定律：对所有 $i$ 求和两边恒等（收入=支出），故 $N$ 个工资方程只有 $N-1$ 个独立，须做归一化（本页代码令 $\sum_i w_i L_i=1$）。
$X_n^j$ 依赖 $TR_n$，$TR_n$ 依赖 $\pi_{ni}^j X_n^j$，$\pi$ 依赖 $c$ 与 $P$，$P$ 又依赖 $c$……这是一个高维不动点。不能"先算完价格再算支出"，必须在给定 $w$ 的内层里把 $\{P,c,\pi,X,TR\}$ 同时迭代收敛，再回到外层更新 $w$。
05 均衡求解：双层不动点迭代
给定 $\{T_i^j,\theta_j,\beta_j,a_{kj},b_j,\tau_{ni}^j,t_{ni}^j,L_i\}$，算法分两层：
与 02-EK 相同：w*=(w_implied/w)^omega 在归一化后会收敛到"残差处处相等"的伪不动点。多部门系统雅可比矩阵更密、谱半径更接近 1，算术阻尼 + 归一化是唯一稳妥做法。
06 Python 完整求解器（已实跑验证）
下面是本页完整求解器。已在本机实跑：4 国 × 3 部门，份额行和=1；跨国关税减半后各国福利上升 1.2%–1.7%，方向与理论一致。
"""
Caliendo & Parro (2015, AER) 多国多部门 EK 量化贸易模型求解器
- N 国 J 部门, 部门间 IO 关联 a[k,j], 部门特定 theta_j, 从价关税返还
- 外层工资不动点, 内层 {P,c,pi,X,TR} 同时迭代
"""
import numpy as np
np.random.seed(20240912) # 固定种子, 结果可复现
def solve_calibrate_parro(N, J, L, T, theta, beta, a, b, tau, t,
tol=1e-9, outer_damp=0.25):
"""T(N,J) theta(J,) beta(J,) a(J,J)[品,部门] b(J,) tau(N,N,J) t(N,N,J)"""
w = np.ones(N)
P, X, TR = np.ones((N, J)), np.ones((N, J)), np.zeros(N)
for outer in range(50000):
for inner in range(2000):
# (CP1) log 单位成本
logc = beta[None, :] * np.log(w)[:, None] \
+ (1 - beta[None, :]) * (np.log(P) @ a)
c = np.exp(logc)
# (CP2) 到岸成本与价格指数
landed = c[None, :, :] * tau * (1.0 + t) # (Ndest,Norig,J)
S = T[None, :, :] * np.power(landed, -theta[None, None, :])
Ssum = S.sum(axis=1)
P_new = Ssum ** (-1.0 / theta)[None, :]
# (CP3) 份额
pi = S / Ssum[:, None, :]
# (CP5) 关税收入 TR[n]=Σ_j X[n,j] Σ_i pi*t/(1+t)
tariff_rate = (pi * t / (1.0 + t)).sum(axis=1)
TR_new = (X * tariff_rate).sum(axis=1)
# (CP4) 可支配收入 + IO 中间品需求
income = w * L + TR_new
inter = (X * (1 - beta)[None, :]) @ a.T # 中间品对j的需求
X_new = b[None, :] * income[:, None] + inter
if np.max(np.abs(P_new / P - 1)) < tol and \
np.max(np.abs(X_new / (X + 1e-300) - 1)) < tol:
P, X, TR = P_new, X_new, TR_new; break
P = 0.5 * P + 0.5 * P_new
X = 0.5 * X + 0.5 * X_new
# (CP6) 工资出清: w_i L_i = Σ_j Σ_n pi*X/(1+t)
sales = (pi / (1.0 + t) * X[:, None, :]).sum(axis=(0, 2))
w_implied = sales / L
w_new = (1 - outer_damp) * w + outer_damp * w_implied
w_new = w_new / np.sum(w_new * L)
if outer > 5 and np.max(np.abs(w_implied / w - 1)) < tol:
w = w_new; break
w = w_new
return dict(w=w, P=P, pi=pi, X=X, TR=TR)
# ---------------- 4 国 x 3 部门示例 ----------------
N, J = 4, 3
L = np.array([1.0, 0.8, 1.2, 0.6])
T = np.array([[1.20,1.00,0.90],[0.90,1.30,0.80],
[0.80,0.90,1.40],[0.70,0.80,1.10]])
theta = np.array([4.0, 5.0, 6.0]) # 部门特定贸易弹性
beta = np.array([0.40, 0.35, 0.30]) # 劳动份额
b = np.array([0.45, 0.35, 0.20]) # 最终消费份额
a = np.array([[0.50,0.25,0.20], # a[k,j]=j用k, 列和=1
[0.30,0.50,0.30],
[0.20,0.25,0.50]])
rng = np.random.RandomState(0)
base = np.array([[1,1.6,1.9,2.3],[1.6,1,1.7,2.0],
[1.9,1.7,1,1.8],[2.3,2.0,1.8,1]])
tau = np.ones((N, N, J)); t = np.zeros((N, N, J))
for j in range(J):
tau[:, :, j] = np.maximum(base * (1 + 0.02*rng.randn(N,N)), 1)
np.fill_diagonal(tau[:, :, j], 1.0)
t[:, :, j] = np.where(np.eye(N, dtype=bool), 0,
rng.uniform(0.05, 0.15, (N, N)))
res = solve_calibrate_parro(N, J, L, T, theta, beta, a, b, tau, t)
print("工资 w =", np.round(res["w"], 4))
print("份额行和(=1) =", round(res["pi"].sum(axis=1).mean(), 6))
# 反事实: 跨国关税减半
t1 = t.copy(); t1[~np.eye(N, dtype=bool)] *= 0.5
res1 = solve_calibrate_parro(N, J, L, T, theta, beta, a, b, tau, t1)
Pidx0 = np.exp((b * np.log(res["P"])).sum(axis=1))
Pidx1 = np.exp((b * np.log(res1["P"])).sum(axis=1))
print("福利变化 % =", np.round(((res1["w"]/Pidx1)/(res["w"]/Pidx0)-1)*100, 3))
基准工资 [0.284 0.297 0.265 0.267]，份额行和 = 1.000000；跨国关税减半后福利变化约 +1.2%~+1.7%。你可以把 $N$ 改到 40、$J$ 改到 30，算法仍然收敛——这正是 CP 用 WIOD（40 国 × 35 部门）的原因。
07 估计与数据：$\theta_j$、IO 表、关税
分部门贸易弹性 $\theta_j$
这是 CP 相对单部门 EK 最关键的估计对象。$\theta_j$ 不能用一个全局值，必须按部门估计。识别来源：
- Caliendo-Parro (2015)：用 1993–2001 年 NAFTA 生效前后关税随时间/部门的变异识别 $\theta_j$，制造业 40 个 ISIC 部门，$\theta_j$ 落在 2.5–10，多数在 4–5。
- 货物 vs 服务：服务部门 $\theta$ 通常很大（≈10+），因为服务"贸易成本"主要是非关税壁垒，价格对成本不敏感。
- Broda & Weinstein (2006)：用行业级替代弹性 $\sigma$ 作 $\theta_j$ 近似（$\theta_j\approx\sigma_j-1$）。
数据来源（具体变量）
| 用途 | 数据库 | 变量 / 口径 |
|---|---|---|
| IO 矩阵 $a_{kj}$、增加值份额 $\beta_j$ | WIOD 2016 / OECD ICIO | 国家×国家×部门中间使用矩阵；列归一得 $a_{kj}$，增加值/总产出得 $\beta_j$。 |
| 部门双边贸易流 $\pi_{ni}^j$ | CEPII BACI / UN Comtrade | HS 6 位 CIF 进口，按 ISIC/BEC 归部门；行归一得份额。 |
| 从价关税 $t_{ni}^j$ | WITS TRAINS / WTO TAO | MFN 与 applied 关税，HS 6 位；CP 用 MFN 构造 $\hat t$。 |
| 冰山成本 $\tau_{ni}^j$ | 结构引力残差（见 10 页） | 估 $\ln X_{ni}^j=\alpha_i+\beta_n-\theta_j\ln\tau$，指数化。 |
| 中国海关微观 | 中国海关数据库 + ASIF | 企业-月度进出口，加工贸易与一般贸易分列，估分部门 $\theta_j$。 |
| 中国 IO 表 | 国家统计局投入产出表 / OECD-ICIO 中国行 | 分 42 部门 IO 系数；做中国多区域版本用多区域 IO（MRIO）。 |
08 福利与帽子代数反事实
多部门 ACR 福利不是单一 $-\theta$ 公式，而是部门加权的：一国 $n$ 的福利（消费等价变化）来自两部分——直接的消费侧，与通过中间品成本下降间接传导的供给侧。
直觉：即使某国自己不生产某进口品，关税下降也能通过"中间品更便宜→本国企业成本下降→工资上升"的链条改善福利。这就是 CP (2015) 发现 NAFTA 福利比单部门 ACR 大出一大截的原因。
帽子代数（详见 05 页）：反事实只需要基准 $\pi_{ni}^{j,0}$、$X_n^{j,0}$、外生 $\hat\tau_{ni}^j$、$\hat t_{ni}^j$，不需要重估 $T_i^j$。份额帽子方程：
本页代码为教学清晰起见，用"绝对空间重解"实现反事实（直接改 $\tau,t$ 再求一遍）。论文里为了速度与数据可得性，改用帽子代数（只需基准份额）。
09 论文案例与常见错误
常见错误（6 条）
$\theta_j$ 必须分部门。用全局 $\theta=5$ 同时套到化工（弹性高）和服务（弹性低），会把关税对制造业的福利效应低估一半。CP 的核心贡献就是分部门 $\theta_j$。
$a_{kj}$ 是"部门 $j$ 用部门 $k$"（列和为 1）。中间品对部门 $j$ 的需求是 $\sum_k a_{jk}(1-\beta_k)X^k$（行=采购品）。写反后，中间品渠道方向颠倒，工资迭代不收敛或福利符号错。本页代码用 (X*(1-beta)) @ a.T 明确"行=采购品"。
$i$ 国真正收到的是 CIF 价值 $\pi_{ni}^jX_n^j/(1+t_{ni}^j)$，不是含关税的 $\pi X$。若直接用 $\pi X$，关税收入被重复计入 $i$ 国收入，工资虚高、关税反事实结论完全错。
家庭可支配收入是 $w_nL_n+TR_n$，不是 $w_nL_n$。漏掉返还，相当于关税被"扔进海"，会系统性高估关税的负福利（实际上关税扭曲只来自相对价格，不来自收入损失）。
$(c\tau(1+t))^{-\theta_j}$ 在 $\theta_j=10$、$\tau=3$ 时约 $3^{-100}\approx10^{-48}$，numpy 下溢成 0，份额变 NaN。必须 np.power 配合 log，或直接在 log 空间做乘法。
若内层 $X,TR$ 还没收敛就跳出去更新 $w$，外层会把"未收敛的 IO 噪声"当成工资信号放大，迭代震荡发散。务必内层 $|X^{new}/X-1|<10^{-9}$ 后才更新外层。
∑ 方程总清单与难度前置
| 编号 | 名称 | 公式 | 节 |
|---|---|---|---|
| (CP1) | 部门单位成本 | $c_i^j=w_i^{\beta_j}(\prod_k P_i^{k\,a_{kj}})^{1-\beta_j}$ | 02 |
| (CP2) | 部门价格指数 | $P_n^j=[\sum_i T_i^j(c_i^j\tau_{ni}^j(1+t_{ni}^j))^{-\theta_j}]^{-1/\theta_j}$ | 03 |
| (CP3) | ★ 部门贸易份额 | $\pi_{ni}^j=T_i^j(c_i^j\tau_{ni}^j(1+t_{ni}^j))^{-\theta_j}/\sum_k\cdots$ | 03 |
| (CP4) | ★ 部门支出 | $X_n^j=b_j(w_nL_n+TR_n)+\sum_k a_{jk}(1-\beta_k)X_n^k$ | 04 |
| (CP5) | 关税收入返还 | $TR_n=\sum_j\sum_i\pi_{ni}^jX_n^j\,t_{ni}^j/(1+t_{ni}^j)$ | 04 |
| (CP6) | ★ 工资出清 | $w_iL_i=\sum_j\sum_n\pi_{ni}^jX_n^j/(1+t_{ni}^j)$ | 04 |
| (CP7) | ★ 多部门福利 | $\hat W_n=\hat w_n/\prod_j(\hat P_n^j)^{b_j}$ | 08 |
| 维度 | 要求 |
|---|---|
| 数学 | Fréchet 闭包（见 02）、$J\times J$ IO 矩阵、双层不动点、log 空间幂次 |
| 计量 | 分部门引力估计 $\theta_j$（见 10 页）、PPML（见 03 页） |
| 编程 | numpy 三维张量 $(N,N,J)$、内外层嵌套迭代；先跑通 02 单部门 EK |
| 前置阅读 | 02-EK → 04-ACR → 05-帽代数 → 本页 12 → 11-GVC |
| 交叉链接 | 02-EK 单部门母版 · 03-引力 PPML · 11-GVC 增加值链条 |