结构引力方程：从 Armington 微观基础到 PPML 估计与 ACR 福利
本页自包含地走完一条完整链条：Anderson-van Wincoop (2003) 的理论推导 → 为什么 OLS 对数化不一致 → PPML 双向固定效应估计 → 贸易弹性识别 → 数据要求 → 边界效应与 RTA 应用 → ACR 福利公式 → 关税反事实。
01 经济环境与假设
结构引力（structural gravity）是国际贸易实证与量化一般均衡之间的桥梁。它不是一个"回归方程"，而是一组从明确偏好、技术、贸易成本假设中推导出的均衡关系。要把它用对，必须先把假设写清楚。
经济由 $N$ 个国家（或地区）组成。核心假设按重要性排列如下：
- 假设 A1：Armington 偏好。消费者偏好差异化产品的"产地"——同样一件衬衫，中国产与越南产在效用上不完全替代。$n$ 国代表性消费者在 CES 聚合函数下消费来自所有来源 $i$ 的品种，替代弹性 $\sigma>1$。这一假设给出了 CES 支出份额这一关键微观基础。
- 假设 A2：完全竞争与 iceberg 贸易成本。每单位产品在 $i$ 国出厂价 $p_i$，运到 $n$ 国必须运输 $\tau_{ni}\geq 1$ 单位才能到达 1 单位，故到岸价 $p_{ni}=\tau_{ni} p_i$。贸易成本 $\tau_{ni}$ 是双边的、对称或不对称均可。
- 假设 A3：贸易平衡。每期 $n$ 国总支出 $X_n$ 等于其总产出 $Y_n$（长期封闭经济下的标准假设；开放模型中可加贸易赤字 $D_n$ 作为外生参数）。
- 假设 A4：同质性偏好与生产。所有国家偏好结构相同（相同 $\sigma$），技术是规模报酬不变。这是把 $\sigma$ 当作"贸易弹性"在跨国间可比的关键。
- 假设 A5：无中间投入（基准版）。结构引力基准模型只含劳动一种要素；加入投入产出表后即变为 Caliendo-Parro (2015) 的多部门版本，引力形式保留但需按部门估计。
1962 年 Tinbergen 把 Newton 引力公式类比到贸易：$X_{ni}=G\,Y_n^a Y_i^b/D_{ni}^\rho$。这一经验方程在数据上拟合极好（$R^2$ 常达 0.6–0.8），但在 1962–1990 年代一直是"有数据无理论"的黑箱。Anderson (1979) 第一次用 Armington 偏好给出微观基础，而 Anderson & van Wincoop (2003, AER) 才真正把它做成"structural gravity"——结构引力，并一举解决了 McCallum (1995) 的 border puzzle。
理解这一页的关键，是把"引力方程"从一个回归技术问题提升到一个结构模型问题：方程里每个系数都对应一个经济参数（$\sigma-1$ 或 $\theta$ 是替代弹性、$\gamma$ 是贸易成本的对数参数），每个残差都有结构解释（多边阻力未被充分控制、未观察的贸易摩擦）。只有这样，估计出的系数才能反过来用于 ACR 福利反事实，而不是停留在"$R^2$ 很高"的描述性统计。
02 完整推导：从 Armington 需求到引力方程
Step 1：CES 需求函数（完整 FOC 推导）
$n$ 国消费者求解：
构造拉格朗日 $\mathcal L=U_n-\lambda(\sum_i p_{ni}c_{ni}-X_n)$。对 $c_{ni}$ 求 FOC：
对两品种 $i,k$ 相除得 $c_{ni}/c_{nk}=(p_{ni}/p_{nk})^{-\sigma}$，代入预算约束并定义价格指数 $P_n^{1-\sigma}=\sum_i p_{ni}^{1-\sigma}$，解出：
$$c_{ni}=\frac{X_n}{P_n}\left(\frac{p_{ni}}{P_n}\right)^{-\sigma}$$一阶条件给出对 $i$ 国产品的支出：
代入到岸价 $p_{ni}=\tau_{ni} p_i$：
Step 2：市场出清与多边阻力定义（双向完整推导）
外向 MR（出口侧）：$i$ 国总产出必须等于全世界对它的购买：$Y_i = \sum_n X_{ni}$。把支出份额代入：
两边除以世界总产出 $Y_w=\sum_n Y_n$，并记世界支出份额 $\theta_n=Y_n/Y_w$：
$$\Pi_i^{1-\sigma}\equiv\sum_n\left(\frac{\tau_{ni}}{P_n}\right)^{1-\sigma}\theta_n,\qquad p_i^{1-\sigma}=\frac{Y_i}{Y_w\Pi_i^{1-\sigma}}=\frac{\theta_i}{\Pi_i^{1-\sigma}}$$内向 MR（进口侧）· 原文"同理可得"逐步写出：把价格指数定义 $P_n^{1-\sigma}=\sum_i(\tau_{ni}p_i)^{1-\sigma}$ 中的 $p_i^{1-\sigma}=\theta_i/\Pi_i^{1-\sigma}$ 代入：
Step 3：结构引力方程（消去出厂价）
把 (G2) 与 $p_i^{1-\sigma}=\theta_i/\Pi_i^{1-\sigma}$ 联立：
代入 $\theta_i=Y_i/Y_w$，整理：
$$\boxed{\;X_{ni} = \dfrac{Y_n\,Y_i}{Y_w}\,\dfrac{\tau_{ni}^{\,1-\sigma}}{\Pi_i^{\,1-\sigma}\,P_n^{\,1-\sigma}}\;}$$这就是 Anderson-van Wincoop (2003) 的核心方程。它与朴素引力的关键区别在于分母上的 $\Pi_i P_n$：双边贸易流不仅取决于双边成本 $\tau_{ni}$，还取决于两国与世界其他国家的平均阻力。
03 多边阻力项 (Multilateral Resistance)
Anderson-van Wincoop (2003) 的核心贡献是认识到：双边贸易流由双边成本 + 多边平均成本共同决定。
- 内向多边阻力 $P_n$：$n$ 国作为进口国，从全世界买东西的平均难易程度；
- 外向多边阻力 $\Pi_i$：$i$ 国作为出口国，把货卖到全世界的平均难易程度。
直觉：加拿大与美国相邻，按理加美贸易应该很多。但若加拿大离所有欧洲国家都很远，加拿大对美国的"相对距离"其实并不远——这会放大加美贸易。McCallum (1995) 曾发现"美国内部贸易比美加跨境贸易大 22 倍"（border puzzle）。Anderson-van Wincoop 证明：一旦控制多边阻力，这个 border effect 从 22 倍缩水到约 10 倍，仍巨大但不再神秘。
$\Pi_n$ 对每个进口国 $n$ 是固定的、$P_i$ 对每个出口国 $i$ 是固定的；在面板中它们随年份变化。因此现代引力回归必须用 进口国×年固定效应 与 出口国×年固定效应 吸收，不能只放 GDP 控制变量。
3.1 MR 方程组的数值迭代求解
结构估计之后若要做一般均衡反事实，需要从估计出的 $\hat\tau_{ni}$ 反解 $\Pi_i, P_n$ 水平值。(G3)(G4) 是一个 $2N$ 维非线性不动点方程组。Anderson-van Wincoop (2003) 建议简单不动点迭代：从 $P_n^{(0)}=1$ 出发，交替更新外向、内向 MR，每步末尾做归一化（MR 只有相对水平有意义，须固定一个标度）：
# ============================================================
# 多边阻力 (MR) 方程组数值求解
# Pi_i^(1-sigma) = sum_n (tau_ni / P_n)^(1-sigma) theta_n
# P_n^(1-sigma) = sum_i (tau_ni / Pi_i)^(1-sigma) theta_i
# 归一化: sum_i theta_i * Pi_i^(1-sigma) = 1 (或固定某国 P=1)
# ============================================================
import numpy as np
def solve_mr(tau, theta, sigma=8.0, tol=1e-10, maxit=10000):
"""
tau: NxN 贸易成本矩阵 tau[n,i] (进口国n从出口国i)
theta: N 维各国支出份额 theta_i = Y_i/Y_w
sigma: 替代弹性
返回: Pi (外向MR), P (内向MR) 两个 N 维向量
"""
N = len(theta)
s = 1.0 - sigma # = -(sigma-1)
Pi = np.ones(N)
P = np.ones(N)
for it in range(maxit):
# 外向 MR: Pi_i^s = sum_n (tau[n,i]/P_n)^s theta_n
Pi_new = np.zeros(N)
for i in range(N):
Pi_new[i] = np.sum(theta * (tau[:, i] / P)**s)
# 内向 MR: P_n^s = sum_i (tau[n,i]/Pi_i)^s theta_i
P_new = np.zeros(N)
for n in range(N):
P_new[n] = np.sum(theta * (tau[n, :] / Pi)**s)
# 归一化 (MR 只确定到一个标度)
scale = np.sum(theta * Pi_new**s)
Pi_new /= scale ** (1.0 / s)
P_new /= scale ** (1.0 / s)
if (np.max(np.abs(Pi_new - Pi)) < tol and
np.max(np.abs(P_new - P)) < tol):
Pi, P = Pi_new, P_new
break
Pi, P = Pi_new, P_new
return Pi, P, it
# ---- 示例: 3 国, 距离型 tau = d_ij^rho ----
N = 3
Y = np.array([1.0, 2.0, 3.0])
theta = Y / Y.sum()
d = np.array([[1.0, 3.0, 5.0],
[3.0, 1.0, 2.0],
[5.0, 2.0, 1.0]])
rho = 0.2
tau = d ** rho
np.fill_diagonal(tau, 1.0)
Pi, P, iters = solve_mr(tau, theta, sigma=8.0)
print(f"迭代收敛于第 {iters} 步")
print("外向 MR Pi:", np.round(Pi, 4))
print("内向 MR P :", np.round(P, 4))
# 贸易量: X_ni = Y_n Y_i / Y_w * tau_ni^(1-sigma) / (Pi_i^(1-sigma) P_n^(1-sigma))
Yw = Y.sum()
X = (Y[:, None] * Y[None, :] / Yw
* tau**(1-sigma)
/ (Pi[None, :]**(1-sigma) * P[:, None]**(1-sigma)))
print("双边贸易份额 X_ni/Y_n 矩阵:")
print(np.round(X / Y[:, None], 3))
MR 方程组对整体标度不变：把 $\Pi_i, P_n$ 同时乘 $c$ 引力方程不变。迭代中每步做 $\sum\theta_i\Pi_i^{1-\sigma}=1$ 归一化即可。若两国对称（相同 $Y$、对称 $\tau$），收敛后 $\Pi_i=P_i$，可用作代码正确性自检。
04 为什么 OLS 错：Jensen 不等式与零贸易流
把 $\ln X_{ni}$ 直接 OLS 看似自然，但 Silva & Tenreyro (2006, REStat) 给出了系统批判：
问题 1：异方差 + Jensen 不等式
真实数据生成过程是水平值形式 $X_{ni} = \exp(\mathbf{x}_{ni}'\beta)\,\eta_{ni}$，其中 $E[\eta_{ni}|\mathbf{x}]=1$。取对数：$\ln X_{ni} = \mathbf{x}_{ni}'\beta + \ln\eta_{ni}$。若 $\eta_{ni}$ 异方差（贸易数据几乎必然如此），则由 Jensen 不等式，$E[\ln\eta_{ni}|\mathbf{x}] = -\tfrac12\,\text{Var}(\ln\eta|\mathbf{x}) \neq 0$，且该偏误依赖于 $\mathbf{x}$。距离越远的观测方差越大，对数变换把其条件均值系统性压低——OLS 距离系数因此被过度负偏（常被估成 $-1.2$ 到 $-1.5$）。
问题 2：零贸易流被丢弃
$\ln 0=-\infty$。在 200×200 国对样本中零贸易流常占 20%–50%。OLS 只能用 $\ln X>0$ 的子样本，造成样本选择偏误。用 $\ln(X+1)$ 或 $\ln(X+c)$ 补救会引入新的、不可解释的偏误。
问题 3：加性误差与乘性误差混淆
引力方程是结构性的乘性关系，OLS 对数化隐含假设误差是加性的、同方差的；这两条都与数据相悖。
05 PPML 估计原理与一致性
Poisson Pseudo-ML 直接在水平值 $X_{ni}$ 上估计条件均值。完整推导：设 $X_{ni}\sim\text{Pois}(\mu_{ni})$，$\mu_{ni}=\exp(\mathbf{x}_{ni}'\beta)$。对数似然：
对 $\beta$ 求导（注意 $\partial\mu/\partial\beta=\mu\cdot\mathbf{x}_{ni}$）：
$$\frac{\partial\ell}{\partial\beta}=\sum_{ni}\left(\frac{X_{ni}}{\mu_{ni}}-1\right)\frac{\partial\mu_{ni}}{\partial\beta} =\sum_{ni}\mathbf{x}_{ni}\left[X_{ni}-\exp(\mathbf{x}_{ni}'\beta)\right]=0$$PPML 一致的原因有三：
- 只依赖一阶矩。PPML 是 pseudo-ML，即使真实分布不是 Poisson，只要条件均值正确设定为 $\exp(\mathbf{x}\beta)$，一阶条件矩就正确，$\hat\beta$ 一致；
- 对零贸易流天然容纳。Poisson 分布在 0 处概率为正，$\Pr(X=0)=\exp(-\lambda)$，零贸易流不再是问题；
- 异方差不破坏一致性。误差结构只影响标准误的估计（用稳健/聚类标准误修正），不影响点估计。
实务上把引力因子写成线性指数形式：$\mathbf{x}_{ni}'\beta = (1-\sigma)\ln\tau_{ni} + \text{MR FE}$，其中 $\ln\tau_{ni} = \rho_D \ln D_{ni} + \gamma_1 \text{contig}_{ni} + \gamma_2 \text{lang}_{ni} + \gamma_3 \text{colony}_{ni} + \gamma_4 \text{RTA}_{ni}$。MR 项由出口国×年、进口国×年固定效应用 dummy 吸收。
06 Stata 代码：ppmlhdfe 双向固定效应
下面是行业事实标准的 Stata 工作流（Silva-Tenreyro 2006；Head-Mayer 2014 手册）。
* ================================================================
* Structural Gravity with PPML —— 标准工作流
* 数据: 双边贸易 flow X_ni、距离、接壤、语言、殖民、RTA
* ================================================================
clear all
set more off
* ssc install ppmlhdfe, replace // 首次运行
use "gravity_data.dta", clear
* 变量: flow(水平值!) exp imp year lndist contig lang colony rta
* ---------- 1. OLS 基准（被批判的对照） ----------
gen ln_flow = log(flow)
reghdfe ln_flow lndist contig lang colony rta, ///
absorb(exp#year imp#year) vce(robust)
estimates store ols_fe
* ---------- 2. PPML 推荐版：ppmlhdfe ----------
* 同时吸收出口-年、进口-年固定效应 = 时变 MR
ppmlhdfe flow lndist contig lang colony rta ///
, absorb(exp#year imp#year) vce(robust)
estimates store ppml_fe
* ---------- 3. 对比系数 ----------
* 距离弹性 PPML 预期 -0.7 ~ -1.0；OLS 常偏到 -1.2 以下
estimates table ols_fe ppml_fe, ///
b(%9.4f) star(.1 .05 .01) keep(lndist contig lang colony rta)
* ---------- 4. 边界效应（含国内贸易 X_ii） ----------
* 需先用 (总产出 - 总出口) 构造国内贸易流
* gen intl = (exp != imp)
* ppmlhdfe flow lndist intl contig lang colony rta, ///
* absorb(exp#year imp#year) vce(robust)
* exp(intl) 显著为负 = 边界效应；预期 exp(intl)~0.1~0.3
* ---------- 5. 聚类标准误 ----------
* 双边数据相关性来自同一 country pair，建议 vce(cluster pair_id)
* ppmlhdfe flow ... , absorb(exp#year imp#year) vce(cluster pair_id)
* ---------- 6. 缺失/分离检查 ----------
* ppmlhdfe 会报告 dropped obs；关注"完美分离"是否合理
(1) ppmlhdfe 比 poisson 快 10–100 倍，对高维固定效应内存友好；(2) 必须 absorb(exp#year imp#year)，否则漏掉时变 MR；(3) 国内贸易流 $X_{ii}$ 用（总产出 − 总出口）构造，不要当 0；(4) 标准误按 pair 聚类，因为同一国家对在多个年份相关。
07 贸易弹性识别与系数解读
引力回归估计的核心结构参数是贸易弹性。在 Armington 设定中它是 $\sigma-1$；在 EK 设定中它是 Fréchet 形状参数 $\theta$。两者在引力方程中扮演完全相同的角色：
识别逻辑：若外部给出 $\theta=5$（Simonovska-Waugh 2014 综述的典型值），距离系数 $-1.0$ 就意味着 $\rho_D=0.2$——距离对贸易成本的弹性。反之若直接用关税变化识别（关税 $\to$ 进口份额），则 $\theta$ 可直接从数据估出。
系数解读要点：
- 距离系数 $\approx -0.8$：距离翻倍，双边贸易下降约 $2^{0.8}-1\approx 74\%$；
- 共同语言 $\approx 0.4$：共同语言使贸易上升 $e^{0.4}-1\approx 49\%$；
- 殖民关系 $\approx 0.6$：前殖民地与宗主国贸易高约 $e^{0.6}-1\approx 82\%$；
- RTA $\approx 0.4$：加入 RTA 使贸易上升约 49%，但该估计对内生性极敏感。
注意：引力方程中 GDP 项被双向固定效应吸收，不应再把 $\ln Y_n, \ln Y_i$ 放进去——否则与固定效应共线。
08 数据要求与变量构造
8.1 双边贸易流
- UN Comtrade：联合国官方双边 HS 码贸易流，免费公开，覆盖 1988 年至今，按 reporter-partner-year 报告；
- CEPII BACI（Gaulier-Zignago 2010）：在 Comtrade 基础上统一报关价、清理重复报告（reporter 与 partner 对同一流量的不一致），是学者事实标准；
- WIOD / OECD-ICIO：世界投入产出表，含国内贸易 $X_{ii}$，做 ACR 反事实必需；
- 国内贸易流构造：$X_{ii} = $ 总产出（GDP 增加值需上调为 gross output，可用投入产出表）− 总出口。
8.2 引力变量（全部来自 CEPII）
- CEPII GeoDist 数据库（Mayer-Zignago 2011）：距离 $D_{ni}$（dist 人口加权距离 vs distcap 首都距离，实证推荐 distw）、接壤 contig、共同语言 lang（ethno 官方语言或语言族）、殖民关系 colony（殖民关系）与 comcol（共同殖民地）；
- Egger-Larch RTA 数据库 / WTO RTA：区域贸易协定生效年份与类型（FTA/关税同盟/经济一体化）；
- 关税：MAcMap-HS6（CEPII/ITC）、TRAINS（World Bank）、WTO IDB，做 $\theta$ 识别用。
8.3 GDP 与宏观控制变量
- Penn World Table (PWT 10.0)：跨国购买力平价 GDP（rgdpna）、就业、资本存量，做跨国长期面板首选；
- IMF World Economic Outlook (WEO)：名义 GDP 与预测值，覆盖最新年份；
- World Bank World Development Indicators (WDI)：GDP、人口、汇率，免费公开，与 Comtrade 对码最方便；
- 注意：结构引力 + 双向固定效应下 GDP 被吸收，不需要在回归里放；但帽子代数反事实需要各国 $Y_i$ 水平值，从上述三库任选其一（口径需一致）。
8.4 中国数据
- 中国海关数据库（2000–至今）：企业-产品(HS8)-目的地(243 国)层面进出口，含交易值、数量、贸易方式（一般/加工），做中国出口引力、企业层面二元边际首选；
- 中国工业企业数据库（1998–2013，规上工业）：与海关库按企业名/邮编+电话匹配（"工企-海关匹配"，如余淼杰、戴觅等用），可同时识别生产率与出口状态；
- 国家统计局 / CEIC：分省 GDP、行业产值、进出口总额宏观序列，用于省级引力（中国省际贸易、一带一路）；
- 国内贸易流 $X_{ii}$ 构造：对中国，用分省投入产出表（国家统计局 / 国务院发展研究中心）或"总产出 − 总出口"，不可当 0。
面板建议：200 个国家 × 30 年 ≈ 60 万观测（其中一半是零）。双向固定效应吸收后，剩余识别来自"随 pair 变化的变量"（距离不随时间变 → 只能靠 cross-section 识别；RTA 生效前后变化靠 within-pair 时间 variation）。
09 应用解读：贸易成本、边界效应、RTA
估计出引力方程后，可以回答三类政策问题：
- 贸易成本估计。把 $\tau_{ni}$ 反推出来：$\hat\tau_{ni} = \exp[(\hat\gamma_1\text{contig}+\cdots)/\hat\theta]$。典型结果：距离每增加 1000 km 约等价于加征 20%–30% 从价税。
- 边界效应。加入"国内贸易 vs 国际贸易"dummy $intl$，其系数 $\exp(\hat\gamma_{intl})$ 即跨境"未观察到的"边界成本。Anderson-van Wincoop 2003 把它从 McCallum 的 22 倍修正到 ~10 倍。
- RTA 效应。加入 RTA dummy，注意内生性：贸易伙伴关系好的国家更可能签 RTA。Baier-Bergstrand (2007) 用 pair FE + IV（历史殖民、外交事件）缓解。
10 福利含义：ACR 公式
PPML 估出的 $\theta$ 是 ACR 公式的"输入参数"。ACR 公式的完整包络定理推导、成立条件、多部门 (Caliendo-Parro) 扩展、常见错误（贸易弹性 vs OLS 距离弹性） → 先学 基础知识库 · 福利分解 (ACR)
结构引力估计的贸易弹性 $\theta$ 可直接接入 Arkolakis-Costinot-Rodriguez-Clare (2012, AER) 福利公式。在贸易平衡 + CES/EK 结构下，$n$ 国从基准到新均衡的福利变化仅取决于两个充分统计量：
其中 $\hat\lambda_{nn} = \lambda'_{nn}/\lambda_{nn}$ 是本国消费中本国产品份额的变化。直觉：贸易开放使本国份额下降（$\hat\lambda_{nn}<1$），实际工资上升。关键结论：无论微观结构是 Armington、EK 还是 Melitz-Pareto，只要满足"垄断竞争 + 单一生产要素 + 贸易平衡"三条假设，福利变化都只取决于 $\theta$ 与本国份额变化——这就是引力估计为什么重要：它给出了 ACR 的"输入参数"。
11 反事实：关税变化的福利效应
考虑一个反事实：某国把从 $i$ 国的关税从 $t$ 降到 $t'$，即 $\hat\tau_{ni}=(1+t')/(1+t)$。用帽子代数（hat algebra）求解：
迭代直到 $\hat w$ 收敛，再由 ACR 算 $\hat W_n=\hat\lambda_{nn}^{-1/\theta}$。
数值示例
假设 $\theta=5$，某大国单边把平均关税从 10% 降到 0%：
- 本国自足份额 $\lambda_{nn}$ 从 0.8 降到 0.75；
- $\hat W_n = (0.75/0.8)^{-1/5} \approx 1.013$——实际工资上升约 1.3%；
- 外国福利变化取决于其贸易份额与多边阻力反馈，需完整求解器。
11.1 完整帽代数求解器（自包含 Python 代码）
把上面的帽子方程写成可运行代码：给定基准份额矩阵 $\pi_{ni}$ 与贸易成本冲击 $\hat\tau_{ni}$，迭代工资不动点 $\hat w$，再由 ACR 输出各国福利变化。下面以"中美互相加征 10% 关税"为例：
# ============================================================
# 结构引力 (Armington/EK/Melitz-Pareto 通用) 关税反事实
# 输入: 基准份额 pi[n,i]=X_ni/Y_n, 贸易成本冲击 tau_hat[n,i]
# 输出: 新份额, 工资变化 w_hat, 各国福利 W_hat
# ============================================================
import numpy as np
def gravity_counterfactual(pi, tau_hat, theta=5.0, tol=1e-9, maxit=2000):
"""
pi: NxN 基准贸易份额 (行归一化: sum_i pi[n,:] = 1)
tau_hat: NxN 贸易成本变化 (1.1 = 成本+10%)
theta: 贸易弹性 = sigma-1 (Armington) 或 k (Melitz-Pareto)
"""
N = pi.shape[0]
# 收入份额 s_i = Y_i / Y_w; 默认等大国
s = np.ones(N) / N
w_hat = np.ones(N) # 工资变化帽, 初值=1
for it in range(maxit):
# Step 1: 给定 w_hat, 每个进口国 n 的分母 D_n
# D_n = sum_k pi[n,k]*tau_hat[n,k]^(-theta)*w_hat[k]^(-theta)
D = np.sum(pi * tau_hat**(-theta) * w_hat[None, :]**(-theta),
axis=1)
# Step 2: 工资不动点 (贸易平衡, 推导见下)
# 新出口额 = sum_n s_n * w_hat_n * pi_hat[n,i]
# = w_hat_i * s_i
# 代入 pi_hat[n,i] = pi[n,i]*tau^-theta*w_hat_i^-theta / D_n,
# 两边乘 w_hat_i^theta:
# w_hat_i^(1+theta) = (1/s_i) sum_n s_n w_hat_n
# * pi[n,i]*tau_hat[n,i]^(-theta) / D_n
new_w = np.ones(N)
for i in range(N):
num = np.sum(s * w_hat * pi[:, i] * tau_hat[:, i]**(-theta)
/ D)
new_w[i] = (num / s[i]) ** (1.0 / (1.0 + theta))
# 归一化: 平均工资不变(只关心相对变化)
new_w /= new_w.mean()
# 阻尼
w_hat_new = 0.5 * w_hat + 0.5 * new_w
if np.max(np.abs(w_hat_new - w_hat)) < tol:
w_hat = w_hat_new
break
w_hat = w_hat_new
A = pi * tau_hat**(-theta) * w_hat[None, :]**(-theta)
pi_hat = A / A.sum(axis=1, keepdims=True)
# ACR: W_hat_i = (新本国份额/基准本国份额)^(-1/theta)
# 注意: 必须用"份额变化帽" diag(pi_hat)/diag(pi), 不是新份额水平
W_hat = (np.diag(pi_hat) / np.diag(pi)) ** (-1.0 / theta)
return pi_hat, w_hat, W_hat
# ---------- 3 国示例: US, CN, ROW; 中美互加 10% 关税 ----------
# pi[n,i]: 行 n 的进口来源结构; 对称份额保证基准贸易平衡
pi = np.array([
[0.70, 0.15, 0.15], # US
[0.15, 0.70, 0.15], # CN
[0.15, 0.15, 0.70], # ROW
])
tau_hat = np.ones((3, 3))
tau_hat[0, 1] = tau_hat[1, 0] = 1.10 # 中美双边贸易成本 +10%
pi_hat, w_hat, W_hat = gravity_counterfactual(pi, tau_hat, theta=5.0)
names = ["US", "CN", "ROW"]
print("国家 福利变化 工资变化 自足份额变化")
for i, nm in enumerate(names):
lam_hat = pi_hat[i, i] / pi[i, i]
print(f"{nm:4s} {W_hat[i]-1:+7.2%} {w_hat[i]-1:+7.2%} {lam_hat:7.3f}")
# 预期: 中美两国福利均下降(贸易成本上升), ROW 因贸易转移略升
单边关税变化会改变贸易条件（terms of trade），大国与小国福利方向相反；ACR 静态公式假设贸易平衡，忽略动态资本积累与劳动力迁移。做大国反事实必须用完整多部门模型（Caliendo-Parro 2015）。上面的工资不动点迭代就是为了让贸易条件内生化——若简单假设 $w_hat=1$（小国开放经济），会系统性高估关税的保护作用。
12 Python 复现
"""
Structural Gravity 的 PPML 估计
- 模拟 N=50 国的双边贸易流
- 用 Poisson GLM 估计引力方程（双向 FE）
"""
import numpy as np
import pandas as pd
import statsmodels.api as sm
np.random.seed(0)
N = 50
countries = [f"C{i}" for i in range(N)]
# ---------- 1. 生成国家层面数据 ----------
gdp = np.exp(np.random.normal(10, 1.5, N))
dist = np.exp(np.random.normal(2, 0.5, size=(N, N)))
np.fill_diagonal(dist, 1.0)
lang = (np.random.uniform(size=(N,N)) < 0.2).astype(int)
colony = (np.random.uniform(size=(N,N)) < 0.05).astype(int)
contig = (np.random.uniform(size=(N,N)) < 0.1).astype(int)
np.fill_diagonal(lang, 1); np.fill_diagonal(colony, 0); np.fill_diagonal(contig, 1)
# ---------- 2. 模拟贸易流（真实系数） ----------
true_beta = dict(dist=-0.9, lang=0.4, colony=0.6, contig=0.3)
rows = []
for i in range(N):
for j in range(N):
eta = (np.log(gdp[i]) + np.log(gdp[j])
+ true_beta["dist"]*np.log(dist[i,j])
+ true_beta["lang"]*lang[i,j]
+ true_beta["colony"]*colony[i,j]
+ true_beta["contig"]*contig[i,j])
flow = np.random.poisson(np.exp(eta - 20))
rows.append(dict(exp=countries[i], imp=countries[j], flow=flow,
lndist=np.log(dist[i,j]), lang=lang[i,j],
colony=colony[i,j], contig=contig[i,j]))
df = pd.DataFrame(rows)
# ---------- 3. OLS（错的，对照） ----------
df["ln_flow"] = np.log(df["flow"] + 1)
ols = sm.ols("ln_flow ~ lndist + lang + colony + contig + C(exp) + C(imp)",
data=df).fit()
print("OLS（含 +1 偏差）距离系数:", round(ols.params["lndist"], 3))
# ---------- 4. PPML（正确） ----------
X = pd.concat([
df[["lndist","lang","colony","contig"]].astype(float),
pd.get_dummies(df["exp"], prefix="exp", drop_first=True).astype(float),
pd.get_dummies(df["imp"], prefix="imp", drop_first=True).astype(float),
], axis=1)
X = sm.add_constant(X)
ppml = sm.GLM(df["flow"], X, family=sm.families.Poisson()).fit()
print("PPML 距离系数:", round(ppml.params["lndist"], 3),
"(真实值 -0.9)")
13 论文案例与常见错误
常见错误（5 条）
异方差下 $\ln X=\mathbf{x}\beta+\ln\eta$ 的 OLS 不一致。距离系数常被估到 $-1.2$ 以下，PPML 后回到 $-0.8$ 到 $-1.0$。
MR 是时变的。必须 absorb(exp#year imp#year)，漏掉时变 $\Pi_i$ 会让距离与 RTA 系数偏估。
丢零贸易流 = 选择偏误；$\ln(X+1)$ 引入不可解释的偏误。PPML 天然容纳 0。
exp#year FE 已经吸收出口国 GDP 时变项，再放 $\ln Y_i$ 完全共线；同理 importer 侧。新手常因复制旧版朴素引力代码而报错。
RTA 是内生选择（关系好的国家更可能签约）。Baier-Bergstrand (2007) 建议用 pair FE + IV，或用 structural gravity 反事实验证。直接 OLS/PPML 的 RTA 系数常被高估 2–3 倍。
进阶资料
- Baier & Bergstrand (2007), "Do FTAs Actually Increase Members' International Trade?", JIE。
- Simonovska & Waugh (2014), "The Trade Elasticity", JIE。
- Behrens et al. (2012), "The Distinct Effects of Trade Policy on Trade", JIE。
∑ 方程总清单 · Equation Summary
本模型共 14 个方程，按推导顺序编号如下。
| 编号 | 名称 | 公式 | 所在节 |
|---|---|---|---|
| (G1) | CES 需求 FOC | $U_n^{1/\sigma}c_{ni}^{-1/\sigma}=\lambda p_{ni}$ | 02 |
| (G2) | Armington 支出份额 | $X_{ni}=(p_{ni}/P_n)^{1-\sigma}X_n$ | 02 |
| (G3) | ★ 外向 MR | $\Pi_i^{1-\sigma}=\sum_n(\tau_{ni}/P_n)^{1-\sigma}\theta_n$ | 02 |
| (G4) | ★ 内向 MR | $P_n^{1-\sigma}=\sum_i(\tau_{ni}/\Pi_i)^{1-\sigma}\theta_i$ | 02 |
| (G5) | ★ 结构引力方程 | $X_{ni}=Y_nY_i/Y_w\cdot\tau_{ni}^{1-\sigma}/(\Pi_i^{1-\sigma}P_n^{1-\sigma})$ | 02 |
| (G6) | 对数线性化 | $\ln X_{ni}=...-\theta\ln\tau_{ni}+\epsilon_{ni}$ | 03 |
| (G7) | 贸易成本展开 | $\ln\tau_{ni}=\rho_D\ln D_{ni}+\rho_{contig}\cdot contig+\rho_{lang}\cdot lang+\rho_{RTA}\cdot RTA$ | 06 |
| (G8) | PPML 似然与 score | $\sum\mathbf{x}_{ni}[X_{ni}-\exp(\mathbf{x}_{ni}'\beta)]=0$ | 05 |
| (G9) | ★ PPML 一阶矩 | $E[X_{ni}\mid\mathbf{x}_{ni}]=\exp(\mathbf{x}_{ni}'\beta)$ | 05 |
| (G10) | 对称双边结构引力 | $X_{ij}/Y_i=(\tau_{ij}/\Pi_iP_j)^{1-\sigma}(Y_j/Y_w)$ | 06 |
| (G11) | 关税反事价格冲击 | $\hat\tau_{ni}=1+t'_{ni}/(1+t_{ni})$ | 08 |
| (G12) | 份额帽子方程 | $\hat\pi_{ni}=(\hat\tau_{ni}\hat w_i)^{1-\sigma}/\sum_k\pi_{nk}(\hat\tau_{nk}\hat w_k)^{1-\sigma}$ | 09 |
| (G13) | 工资不动点 | $\hat w_i=\sum_n\pi_{ni}\hat\pi_{ni}\hat w_n w_nL_n/(w_iL_i)$ | 09 |
| (G14) | ★ ACR 福利公式 | $\hat W_n=\hat\pi_{nn}^{-1/(\sigma-1)}$ | 07 |