Eaton-Kortum 模型：Fréchet 生产率分布与引力结构

01 模型设定：Ricardian 世界与连续品

Eaton & Kortum (2002, Econometrica) 把 Ricardo 的比较优势思想做成了可量化模型。世界有 $N$ 个国家。消费者消费一单位连续品 $j \in [0,1]$ 的 CES 篮子：

不同国家生产同一品种 $j$ 的技术不同：$i$ 国在品种 $j$ 上的生产率记为 $z_i(j)>0$（单位劳动产出 $z_i(j)$ 单位产品）。$i$ 国工资为 $w_i$，则生产 $j$ 品种的单位成本为 $w_i/z_i(j)$。iceberg 贸易成本 $\tau_{ni}\geq 1$ 把 $i$ 国产品运到 $n$ 国的到岸成本放大 $\tau_{ni}$ 倍：

完全竞争 + 一价定律：$n$ 国消费者在所有可能的来源地中选最便宜的。于是 $n$ 国实际购买 $j$ 的价格为：

关键问题：$z_i(j)$ 是一个随机变量，由 Fréchet 分布刻画。我们要在这种不确定性下，算出"从 $i$ 国买 $j$"的概率，再对 $j$ 积分，就得到贸易份额。

02 Fréchet 生产率分布

EK 假设各国生产率 $z_i(j)$ 是独立同分布的随机变量，服从两参数 Fréchet 分布：

两个参数各有直观含义：

• $T_i$ 绝对生产率：$T_i$ 越大，分布整体右移——$i$ 国在大多数品种上生产率更高。它对应 Ricardo 模型里的"绝对优势"。
• $\theta$ 比较优势离散度：$\theta$ 越小，分布尾部越厚，各国生产率差异越大——这是"比较优势的强度"。$\theta$ 同时是后面的贸易弹性（trade elasticity）：贸易成本下降 1%，贸易流对它的弹性就是 $\theta$。

密度函数为 $f_i(z) = F_i'(z) = \theta T_i z^{-\theta-1} e^{-T_i z^{-\theta}}$。这个分布有个特殊的数学性质：独立 Fréchet 变量的最小值（按比例缩放后）仍然是 Fréchet。这是整个 EK 推导能闭合的关键。

03 从分布到价格：第一步推导

我们先把价格也变成随机变量。定义 $i$ 国供给 $n$ 国的"单位成本"为：

由于 $z_i(j)$ 是 Fréchet，我们可以做变量替换求 $x_{ni}(j)$ 的分布。记 $c_i \equiv w_i \tau_{ni}$，则 $z_i(j) = c_i / x$。于是：

令 $S_{ni} \equiv T_i (w_i \tau_{ni})^{-\theta}$（一个与品种 $j$ 无关的常数，只依赖国家对 $(n,i)$），则上式简化为：

注意：$x_{ni}(j)$ 本身在 $j$ 上是独立同分布的——所有品种 $j$ 的价格分布都由 $S_{ni}$ 完全决定。这就是为什么模型能"聚合"。

04 贸易份额公式：完整推导

现在求"$n$ 国从 $i$ 国购买品种 $j$"的概率 $\pi_{ni}$。这要求两件事同时成立：(1) $i$ 国的供给价格 $x_{ni}(j)$ 等于某 $p$；(2) 其他所有国家 $k\neq i$ 的供给价格都严格大于 $p$。

Step 1：写出概率密度与生存函数

对任意价格水平 $p$，$i$ 国报价落在 $[p, p+dp]$ 的概率为 $dG_{ni}(p) = S_{ni}\theta p^{\theta-1} e^{-S_{ni} p^\theta} dp$。

其他国家 $k\neq i$ 报价 $>p$ 的概率为 $\Pr(x_{nk}>p) = 1 - G_{nk}(p) = e^{-S_{nk} p^\theta}$。各国独立，所以"所有 $k\neq i$ 都比 $p$ 贵"的概率为：

Step 2：积分

把两件事乘起来，对 $p$ 积分：

令 $\bar{S}_n \equiv \sum_{k=1}^N S_{nk} = \sum_k T_k(w_k\tau_{nk})^{-\theta}$，做变量替换 $u = \bar{S}_n p^\theta$，则 $du = \bar{S}_n \theta p^{\theta-1} dp$，即 $\theta p^{\theta-1} dp = du/\bar{S}_n$。积分变为：

把 $S_{ni}$ 和 $\bar{S}_n$ 代回，得到教科书公式：

这是 QSGE 里最重要的公式之一。注意它与 Armington 模型在结构上几乎一模一样：把 $T_i w_i^{-\theta}$ 类比为"国家 $i$ 的质量"，把 $\tau_{ni}^{-\theta}$ 类比为"贸易成本的引力因子"，分子就是"从 $i$ 国买的吸引力"，分母是全世界的"总吸引力"。但关键差别在于：

• Armington 里弹性是 $1-\sigma$（来自偏好），EK 里弹性是 $-\theta$（来自技术分布）；
• Armington 中每国只生产一种产品，EK 中每国在所有品种上竞争，只是"赢"的概率随 $T_i$ 与 $\tau_{ni}$ 变化；
• EK 自然解释了为什么"贸易同时发生在很多行业"——因为每个品种上的赢家可能不同。

05 价格指数与工资方程

我们还可以推 $n$ 国的 CES 价格指数。在连续品框架下，最便宜价格 $p_n(j)=\min_i x_{ni}(j)$ 的分布函数为：

这仍然是 Fréchet。下面完整推导 CES 价格指数 $P_n$（不引用"标准结果"）。

价格指数的完整推导

$n$ 国消费者在连续品 $j\in[0,1]$ 上的 CES 价格指数为：

由于 $j$ 连续且独立同分布，积分等于对随机变量 $p_n(j)$ 的期望：$\int_0^1 p_n(j)^{1-\sigma}dj=\mathbb{E}[p_n^{1-\sigma}]$。已知 $p_n$ 的分布为 $G_n(p)=1-e^{-\bar S_n p^\theta}$，密度 $g_n(p)=\bar S_n\theta p^{\theta-1}e^{-\bar S_n p^\theta}$。于是：

把 $\pi_{nn}$ 与 $P_n$ 联立，逐步推出用国内份额表示的价格指数：

实际工资为：

与 Armington 完全同构——只是把 $\sigma-1$ 换成了 $\theta$。这正是 ACR (2012) 跨模型同构结论的核心。

从水平公式到 ACR 反事实福利公式

定义帽子 $\hat{x}=x'/x$。由 $P_n=\gamma w_n(T_n/\pi_{nn})^{1/\theta}$（$\tau_{nn}=1$），两边取帽子：

若反事实中技术不变（$\hat{T}_n=1$），则福利变化为：

ACR 公式在 EK 下成立的假设

1. Fréchet 生产率分布（形状参数 $\theta$ 跨国家、跨品种常数）；
2. 贸易平衡：$w_i L_i = \sum_n \pi_{ni} w_n L_n$，无外生贸易赤字；
3. 单一生产要素（劳动）+ 完全竞争 + 规模报酬不变；
4. 反事实中 $T_n$ 不变（只冲击 $\tau$ 或关税）。若 $\hat{T}_n\neq 1$，福利公式要加一项 $\hat{T}_n^{1/\theta}$。

工资决定方程（贸易平衡 · 完整写出）

06 与 Armington 的对比

07 反事实求解：EK 帽子代数

与 Armington 一样，EK 反事实也用精确帽子代数（Dekle-Eaton-Kortum 2007）。但 EK 多了一个外生技术参数 $T_i$，所以帽子方程里多一项 $\hat{T}_i$。

第一步：份额的帽子方程（逐步推导）

基准份额：$\pi_{ni}=T_i(w_i\tau_{ni})^{-\theta}/\sum_k T_k(w_k\tau_{nk})^{-\theta}$。反事实份额：$\pi'_{ni}=T'_i(w'_i\tau'_{ni})^{-\theta}/\sum_k T'_k(w'_k\tau'_{nk})^{-\theta}$。两式相除：

若反事实只冲击贸易成本/关税（$\hat T_i=1$），简化为：

第二步：贸易平衡的帽子方程（逐步推导）

基准 $w_iL_i=\sum_n\pi_{ni}w_nL_n$。反事实 $\hat w_iw_iL_i=\sum_n\pi'_{ni}\hat w_nw_nL_n$（$L_i$ 不变）。把 $\pi'_{ni}=\pi_{ni}\hat\pi_{ni}$ 代入并除以基准 $w_iL_i$：

注意分母是 $w_i L_i$（基准支出），分子里 $\pi_{ni}$ 是基准份额（数据），$\hat{\pi}_{ni}$ 由第一步算出，$\hat{w}_n$ 是当前迭代的工资帽子。

第三步：Jacobi 阻尼迭代算法

第四步：福利与反事实份额

收敛后 $\pi_{ni}'=\pi_{ni}\hat{\pi}_{ni}$，ACR 福利公式直接给出：

若反事实含技术冲击 $\hat{T}_n$，则福利公式修正为 $\hat{W}_n=\hat{T}_n^{1/\theta}\cdot\hat{\pi}_{nn}^{-1/\theta}$（技术进步通过 $\hat{T}_n$ 项直接降低价格）。

08 Python 完整求解与反事实代码

"""
Eaton-Kortum 多国模型：基准求解 + 精确帽子代数反事实
- Part A: 给定 T, L, tau 求解基准均衡 w, P, pi
- Part B: 精确帽子代数反事实（只需基准 pi + 外生 tau_hat）
- Part C: 验证 ACR 福利公式 W_hat = pi_nn_hat^(-1/theta)
"""
import numpy as np

# ---------- 1. 基准均衡求解（绝对空间） ----------
def solve_ek_abs(T, L, tau, theta, tol=1e-10, max_iter=50000, damp=0.3):
    N = len(L)
    w = np.ones(N)
    for it in range(max_iter):
        cost = w[None, :] * tau
        S = T[None, :] * (cost ** (-theta))
        pi = S / S.sum(axis=1, keepdims=True)
        sales = pi.T @ (w * L)
        excess = sales / (w * L)
        w_new = w * excess ** damp
        w_new = w_new / np.sum(w_new * L)
        if np.max(np.abs(excess - 1.0)) < tol:
            w = w_new; break
        w = w_new
    cost = w[None, :] * tau
    S = T[None, :] * (cost ** (-theta))
    pi = S / S.sum(axis=1, keepdims=True)
    P = S.sum(axis=1) ** (-1.0 / theta)   # 忽略常数 gamma
    return w, P, pi

# ---------- 2. EK 精确帽子代数反事实 ----------
def solve_ek_hat(pi0, X0, tau_hat, T_hat=None, theta=8.28,
                tol=1e-10, max_iter=100000, damp=0.2):
    """
    输入:
      pi0     (N,N) 基准贸易份额
      X0      (N,)  基准总支出 w_n L_n
      tau_hat (N,N) 贸易成本冲击的帽子
      T_hat   (N,)  技术冲击帽子，默认全 1
    输出:
      w_hat, pi_hat, pi1 (反事实份额)
    注意: 用算术阻尼 w_new=(1-damp)*w+damp*impl；
          几何阻尼 w*(impl/w)^damp 会收敛到伪不动点（见正文常见错误）。
    """
    N = len(X0)
    if T_hat is None:
        T_hat = np.ones(N)
    w_hat = np.ones(N)
    for it in range(max_iter):
        # 份额帽子方程
        num = T_hat[None, :] * (w_hat[None, :] * tau_hat) ** (-theta)
        den = (pi0 * num).sum(axis=1, keepdims=True)
        pi_hat = num / den
        # 贸易平衡帽子方程
        sales1 = (pi0 * pi_hat).T @ (w_hat * X0)
        w_hat_implied = sales1 / X0
        # 算术阻尼更新（关键：不要用几何阻尼）
        w_hat_new = (1.0 - damp) * w_hat + damp * w_hat_implied
        w_hat_new = w_hat_new / np.sum(w_hat_new * X0)
        if np.max(np.abs(w_hat_implied / w_hat - 1.0)) < tol:
            w_hat = w_hat_new; break
        w_hat = w_hat_new
    pi1 = pi0 * pi_hat
    return w_hat, pi_hat, pi1

# ---------- 3. 数据 ----------
N = 3
T = np.array([1.0, 0.8, 0.6])
L = np.array([1.0, 1.0, 1.0])
tau = np.array([
    [1.0, 1.5, 2.0],
    [1.5, 1.0, 1.8],
    [2.0, 1.8, 1.0],
])
theta = 8.28   # EK (2002) 经典估计

# ---------- 4. 基准 ----------
w, P, pi = solve_ek_abs(T, L, tau, theta)
X = w * L
print("=" * 55)
print("基准 EK 均衡")
print(f"  工资 w          = {np.round(w, 4)}")
print(f"  实际工资 w/P     = {np.round(w/P, 4)}")
print(f"  国内份额 pi_nn   = {np.round(np.diag(pi), 3)}")

# ---------- 5. 反事实：跨国贸易成本下降 10%（帽子代数） ----------
tau_hat = np.ones((N, N))
tau_hat[~np.eye(N, dtype=bool)] = 0.9   # 跨国 tau 下降 10%
w_hat, pi_hat, pi1 = solve_ek_hat(pi, X, tau_hat, theta=theta)

# ACR 福利公式
W_hat_ACR = (np.diag(pi1) / np.diag(pi)) ** (-1.0 / theta)

print("\n" + "=" * 55)
print("反事实：跨国贸易成本下降 10%（帽子代数）")
print(f"  工资帽子 w_hat        = {np.round(w_hat, 4)}")
print(f"  国内份额: {np.round(np.diag(pi),3)} -> {np.round(np.diag(pi1),3)}")
print(f"  ACR 福利变化 W_hat    = {np.round(W_hat_ACR, 4)}")
print(f"  对应福利变化率         = {np.round((W_hat_ACR-1)*100, 3)} %")

# 对照：绝对空间重解
tau_cf = tau.copy()
tau_cf[~np.eye(N, dtype=bool)] *= 0.9
w2, P2, pi2 = solve_ek_abs(T, L, tau_cf, theta)
W_check = (w2/P2) / (w/P)
print(f"  绝对空间重解 W 变化率  = {np.round((W_check-1)*100, 3)} %  (应与 ACR 一致)")

# ---------- 6. 技术赶超反事实：国2 (Foreign1) T 上升 20% ----------
tau_hat2 = np.ones((N, N))
T_hat2 = np.ones(N); T_hat2[1] = 1.2
w_hat2, pi_hat2, pi1b = solve_ek_hat(pi, X, tau_hat2, T_hat=T_hat2, theta=theta)
# 含技术冲击的福利公式: W_hat = T_n_hat^(1/theta) * pi_nn_hat^(-1/theta)
W_hat_tech = T_hat2 ** (1.0/theta) * (np.diag(pi1b)/np.diag(pi)) ** (-1.0/theta)
print("\n" + "=" * 55)
print("反事实：国2 技术 T_hat=1.2")
print(f"  ACR 福利变化 W_hat    = {np.round(W_hat_tech, 4)}")