Melitz (2003) 异质企业贸易模型：从假设到福利与反事实

# ============================================================
# Melitz (2003) 单国模型求解器（Pareto 分布）
# 步骤: ZCP + FE 联立 -> phi* -> phi~ -> M, M_e -> 福利 W
# 反事实: 改变 f_e (贸易自由化等效) 重新求解并量化选择效应
# ============================================================
import numpy as np
from scipy.optimize import brentq

def solve_melitz(sigma=3.8, k=3.4, phi_min=0.2,
                 f=0.1, f_e=1.0, delta=0.02, L=1.0):
    """
    单国 Melitz (Pareto) 均衡求解。
    返回 phi*, phi~, M, M_e, pi_bar, W(福利=实际工资)。
    关键恒等式:
      ZCP 平均利润(常数) = f / (k/(k-(sigma-1)) - 1) - f
      FE  平均利润      = delta * f_e / (1-G(phi*))
      W = M^(1/(sigma-1)) * phi~
    """
    assert k > sigma - 1, "需要 k > sigma-1 才有有限平均生产率"

    # 1-G(phi*) 在 Pareto 下的闭式
    def surv(phi_star):
        return (phi_min / phi_star) ** k if phi_star >= phi_min else 1.0

    # ZCP 给出的平均利润(Pareto 下为常数)
    ratio = k / (k - (sigma - 1))           # (phi~/phi*)^(sigma-1)
    pi_bar_zcp = f / (ratio - 1.0) - f

    # FE: 平均利润 = delta*f_e / surv(phi*)
    # 令 pi_bar_zcp = delta*f_e / surv(phi*), 解 phi*
    # surv(phi*) = (phi_min/phi*)^k = delta*f_e / pi_bar_zcp
    surv_eq = delta * f_e / pi_bar_zcp
    phi_star = phi_min * (surv_eq) ** (-1.0 / k)

    # 加权平均生产率
    phi_tilde = phi_star * (ratio) ** (1.0 / (sigma - 1))

    # 存活企业数 M 与进入者 M_e
    M = L / (sigma * f * ratio)
    M_e = M * delta / surv(phi_star)

    # 福利 W = M^(1/(sigma-1)) * phi~  (w=1)
    W = M ** (1.0 / (sigma - 1)) * phi_tilde

    return dict(phi_star=phi_star, phi_tilde=phi_tilde,
                M=M, M_e=M_e, pi_bar=pi_bar_zcp,
                surv=surv(phi_star), W=W)


def counterfactual(base, **new_params):
    """在 base 基础上改参数，返回福利变化率与 cutoff 变化。"""
    p = dict(sigma=3.8, k=3.4, phi_min=0.2, f=0.1,
             f_e=1.0, delta=0.02, L=1.0)
    p.update(new_params)
    cf = solve_melitz(**p)
    return {
        "phi_star_change": cf["phi_star"] / base["phi_star"] - 1,
        "phi_tilde_change": cf["phi_tilde"] / base["phi_tilde"] - 1,
        "M_change": cf["M"] / base["M"] - 1,
        "welfare_change": cf["W"] / base["W"] - 1,
        "surv_change": cf["surv"] / base["surv"] - 1,
    }


if __name__ == "__main__":
    base = solve_melitz()
    print("=== 基准均衡 ===")
    for k, v in base.items():
        print(f"  {k:<10s} = {v:.4f}")

    # 反事实1: 贸易自由化 -> 等效进入成本 f_e 下降 20%
    cf1 = counterfactual(base, f_e=0.8)
    print("\n=== 反事实: f_e 下降 20% (贸易自由化) ===")
    for k, v in cf1.items():
        print(f"  {k:<16s} = {v:+.2%}")

    # 反事实2: 外生退出率 delta 下降 (更稳定的商业环境)
    cf2 = counterfactual(base, delta=0.01)
    print("\n=== 反事实: delta 0.02 -> 0.01 ===")
    for k, v in cf2.items():
        print(f"  {k:<16s} = {v:+.2%}")

# ============================================================
# 通用分布（对数正态）Melitz 数值求解器
# 残差 H(phi*) = ZCP平均利润(phi*) - FE平均利润(phi*)
# ============================================================
import numpy as np
from scipy.optimize import brentq
from scipy.stats import lognorm
from scipy.integrate import quad

def solve_melitz_general(sigma=3.8, s=0.5, mu=0.0,
                         f=0.1, f_e=1.0, delta=0.02, L=1.0):
    """
    对数正态生产率分布: ln(phi) ~ N(mu, s^2)
    不依赖 Pareto 闭式解, 用 brentq 数值求 ZCP=FE。
    """
    # 对数正态分布 g(phi), CDF G(phi)
    def g(phi):  return lognorm.pdf(phi, s, scale=np.exp(mu))
    def G(phi):  return lognorm.cdf(phi, s, scale=np.exp(mu))

    def tilde_ratio(phi_star):
        """(phi~/phi*)^(sigma-1) = E[phi^(sigma-1)|phi>=phi*] / phi*^(sigma-1)"""
        numer, _ = quad(lambda x: x**(sigma-1) * g(x), phi_star, np.inf)
        return numer / (1.0 - G(phi_star)) / phi_star**(sigma-1)

    def H(phi_star):
        # ZCP 平均利润 (w=1)
        pi_zcp = f * (tilde_ratio(phi_star) - 1.0)
        # FE 平均利润
        pi_fe  = delta * f_e / max(1.0 - G(phi_star), 1e-12)
        return pi_zcp - pi_fe

    # 搜索区间: 分布的 [phi_1%, phi_99%]
    lo = lognorm.ppf(0.01, s, scale=np.exp(mu))
    hi = lognorm.ppf(0.999, s, scale=np.exp(mu))
    phi_star = brentq(H, lo*1.01, hi*0.99, xtol=1e-8)
    # 加权平均生产率
    tilde = (tilde_ratio(phi_star) * phi_star**(sigma-1))**(1/(sigma-1))
    M = L / (sigma * f * tilde_ratio(phi_star))
    W = M ** (1/(sigma-1)) * tilde
    return phi_star, tilde, M, W

phi_star, tilde, M, W = solve_melitz_general()
print(f"对数正态下: phi*={phi_star:.4f}, phi~={tilde:.4f}, M={M:.4f}, W={W:.4f}")
# 与 Pareto 闭式解对比: 同 sigma 下 cutoff 位置接近但尾部行为不同

# ============================================================
# 多国 Melitz-Pareto / EK 统一帽代数求解器
# 方程:
#   pi_hat_ni = tau_hat_ni^(-theta) * (w_hat_i)^(-theta)
#               * pi_ni / sum_k[ tau_hat_nk^(-theta) w_hat_k^(-theta) pi_nk ]
#   w_hat_i 的不动点: sum_n Y_n_hat pi_hat_ni = w_hat_i L_i (贸易平衡)
# 福利: W_hat_i = pi_hat_ii^(-1/theta)
# ============================================================
import numpy as np

def hat_algebra(pi, tau_hat, theta=5.0, s=None, tol=1e-9, maxit=5000):
    """
    pi: NxN 基准贸易份额矩阵, pi[n,i] = 进口国n从i国购买占其总支出份额
    tau_hat: NxN 贸易成本变化 (1=不变, 1.1=成本上升10%)
    theta: 贸易弹性 = Melitz-Pareto k
    s: N 维各国收入份额 Y_i/Y_w; 默认等大国
    返回: pi_hat, w_hat, W_hat(各国福利变化)
    """
    N = pi.shape[0]
    if s is None:
        s = np.ones(N) / N
    w_hat = np.ones(N)

    for it in range(maxit):
        # Step 1: 每个进口国 n 的分母 D_n
        D = np.sum(pi * tau_hat**(-theta)
                   * w_hat[np.newaxis, :]**(-theta), axis=1)

        # Step 2: 工资不动点 (贸易平衡 + 份额方程联立):
        # 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[np.newaxis, :]**(-theta)
    pi_hat = A / A.sum(axis=1, keepdims=True)
    # ACR: W_hat_i = (新自足份额/基准自足份额)^(-1/theta)
    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 = np.array([
    # rows=进口国, cols=出口国:   US    CN   ROW
    [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 = hat_algebra(pi, tau_hat, theta=5.0)
names = ["US", "CN", "ROW"]
print("国家   福利变化  工资变化")
for i, nm in enumerate(names):
    print(f"{nm:4s}  {W_hat[i]-1:+.2%}   {w_hat[i]-1:+.2%}")
# 直觉: 中美两国福利均下降(真实贸易成本上升), ROW 因贸易转移略升