全球价值链与增加值贸易：投入产出、Leontief 逆与 KWW 分解

# ============================================================
# GVC / 增加值贸易核算 (Leontief 逆 + VAX + KWW 分解)
# 参考: Johnson-Noguera (2012); Koopman-Wang-Wei (2014, AER)
# 仅依赖 numpy; 替换 Z,F 为真实 WIOD 即可
# ============================================================
import numpy as np

def build_wiot(G=3, S=2, seed=1):
    """构造模拟世界投入产出表: G国×S部门"""
    rng = np.random.default_rng(seed)
    n = G * S
    # Z[i,j] = 中间品从单元 i 流到单元 j (j 生产用)
    Z = rng.uniform(0.1, 1.0, (n, n))
    # F[i,g] = 最终品从单元 i 流到国家 g
    F = rng.uniform(0.5, 2.0, (n, G))
    x = Z.sum(axis=1) + F.sum(axis=1)          # 总产出(行和)
    A = Z / x[None, :]                          # 直接消耗系数 A[i,j]=Z[i,j]/x_j
    Lmat = np.linalg.inv(np.eye(n) - A)         # Leontief 逆 (I-A)^-1
    va = x - Z.sum(axis=0)                      # 总增加值(列: 产出-中间投入)
    v = va / x                                  # 增加值率 v[j]=VA_j/x_j
    return Z, F, x, A, Lmat, v

def gvc_accounting(Z, F, G, S):
    """核心核算: VAX矩阵 + 总出口 + DVA/FVA + 参与度"""
    n = G * S
    x = Z.sum(axis=1) + F.sum(axis=1)
    A = Z / x[None, :]
    Lmat = np.linalg.inv(np.eye(n) - A)
    va = x - Z.sum(axis=0); v = va / x
    idx = np.array([[g, s] for g in range(G) for s in range(S)])  # 每个单元所属国家

    def mask(g):
        return idx[:, 0] == g                    # 单元 g 的行索引

    # --- 1. VAX[g,m] = g国增加值被m国最终需求吸收 ---
    VAX = np.zeros((G, G))
    for g in range(G):
        for m in range(G):
            out = Lmat @ F[:, m]                 # 满足m最终需求所需总产出
            VAX[g, m] = v[mask(g)] @ out[mask(g)]

    # --- 2. 总出口 EX[g]: 卖给国外的中间品 + 最终品 ---
    EX = np.zeros(G)
    for g in range(G):
        rows = mask(g)
        foreign_col = ~np.isin(idx[:, 0], [g])  # 国外部门
        mid = Z[np.ix_(rows, foreign_col)].sum()
        fin = F[rows][:, [m for m in range(G) if m != g]].sum()
        EX[g] = mid + fin

    # --- 3. g国按来源部门的出口向量 e[g,i] ---
    e = np.zeros((G, n))
    for g in range(G):
        for i in range(n):
            if idx[i, 0] == g:
                mid = Z[i, ~np.isin(idx[:, 0], [g])].sum()
                fin = F[i, [m for m in range(G) if m != g]].sum()
                e[g, i] = mid + fin

    # --- 4. KWW 分解: DVA / FVA / 重复计算 ---
    DVA = np.zeros(G); FVA = np.zeros(G)
    for g in range(G):
        b = Lmat @ e[g]                        # 为生产g出口所需总产出
        for h in range(G):
            va_h = v[mask(h)] @ b[mask(h)]
            if h == g:
                DVA[g] = va_h                   # 本国增加值
            else:
                FVA[g] += va_h                  # 外国增加值
    DC = EX - DVA - FVA                        # 重复计算(会计恒等式)

    # --- 5. GVC 参与度 ---
    vax_ratio = np.array([(VAX[g, :].sum() - VAX[g, g]) / EX[g] for g in range(G)])
    backward = FVA / EX                        # 后向: 外国增加值/总出口
    forward = np.array([(VAX[g, :].sum() - VAX[g, g] - DVA[g] + FVA[g]) / EX[g]
                        for g in range(G)])  # 前向(间接增加值出口)近似
    return dict(VAX=VAX, EX=EX, DVA=DVA, FVA=FVA, DC=DC,
                vax_ratio=vax_ratio, backward=backward, forward=forward)

# ------------------------------------------------------------
if __name__ == "__main__":
    G, S = 3, 2
    Z, F, x, A, Lmat, v = build_wiot(G, S)
    res = gvc_accounting(Z, F, G, S)
    print("VAX 矩阵 (行g=来源国, 列m=吸收国):")
    print(np.round(res["VAX"], 2))
    print("\n总出口 EX:", np.round(res["EX"], 2))
    print("DVA 国内增加值:", np.round(res["DVA"], 2))
    print("FVA 国外增加值:", np.round(res["FVA"], 2))
    print("DC 重复计算:", np.round(res["DC"], 4))
    print("校验 EX = DVA+FVA+DC:", np.allclose(res["EX"], res["DVA"] + res["FVA"] + res["DC"]))
    print("\n--- GVC 参与度 ---")
    for g in range(G):
        print(f"国家{g}: VAX率={res['vax_ratio'][g]:.3f} "
              f"后向参与={res['backward'][g]:.3f}")