BLP 随机系数模型：从 Logit 到 Berry-Levinsohn-Pakes

01 为什么需要随机系数：Logit 的 IIA 问题

上一页我们讲了标准的 logit 离散选择模型。它形式简单、似然函数解析可导，但在差异化产品市场（differentiated product market）的需求估计中存在一个致命缺陷：IIA 性质（Independence of Irrelevant Alternatives，无关替代品独立性）。

回忆一下标准 logit 的选择概率：

其中 $j=0$ 代表外部产品（outside option，即"不购买任何产品"）。IIA 性质意味着：任意两个产品 $j$ 与 $k$ 的相对选择概率之比，只取决于它们自己的效用，与其他产品无关：

著名的"红巴士/蓝巴士"（Red Bus / Blue Bus）悖论

McFadden 当年用一个思想实验批评 IIA：假设一个人在"开车"和"坐红巴士"之间选择，概率各为 1/2。现在加入"蓝巴士"——它和红巴士除了颜色之外完全一样。按 IIA，选择概率必须满足"红巴士:蓝巴士 = 1:1"，而"开车"这一份额必须不变，于是三者各占 1/3。但直觉上，理性人会把原来给红巴士的一半分给蓝巴士，"开车"仍应占 1/2，红/蓝巴士各占 1/4。IIA 强制了"比例替代"（proportionate substitution），这在差异化产品市场里完全不合理：价格上升的产品，其失去的份额应该按"产品之间的相似程度"再分配，而不是机械地按当前份额比例扩散到所有其他产品。

在汽车市场上的具体后果

考虑 Berry, Levinsohn, Pakes (1995) 研究的汽车市场。如果用标准 logit 估计需求，那么当某款"大众甲壳虫"涨价 10%，它流失的份额会按所有其他车（包括保时捷 911、本田思域、雪佛兰索罗德皮卡）的当前份额比例分摊。这显然荒谬：甲壳虫的用户更可能转去 Mini Cooper 或 Fiat 500，而不是转去一辆重型皮卡。标准 logit 的 own-price elasticity（自身价格弹性）只取决于市场份额本身，cross-price elasticity（交叉价格弹性）被 IIA 机械地决定，完全无法反映产品在特征空间中的邻近性。

解决思路很直接：允许消费者对产品特征的偏好是异质的（random coefficients）。这样，喜欢"小尺寸、低油耗"的消费者在甲壳虫涨价时会集中流向类似特征的 Mini，而不是均匀分散。这就是 BLP（1995）的核心思想。

02 BLP 模型设定与效用函数

2.1 消费者效用的随机系数形式

考虑市场 $t$ 中有 $J_t$ 个差异化产品，市场中有 $I$ 个潜在消费者（aggregated 后对应总市场规模 $M_t$）。消费者 $i$ 对产品 $j$ 在时期 $t$ 的间接效用写为：

其中：

• $x_{jt}$：产品 $j$ 的观察特征向量（如马力、油耗、尺寸）；
• $p_{jt}$：价格（内生，因为企业定价时会考虑 $\xi_{jt}$）；
• $\xi_{jt}$：产品 $j$ 的未观察特征（unobserved product characteristic），是结构性误差项，企业和消费者都能看到，但计量经济学家看不到；
• $\varepsilon_{ijt}$：i.i.d. 的 Type-I Extreme Value 误差（和标准 logit 一样）；
• $\beta_i, \alpha_i$：消费者层面的随机系数，随 $i$ 变化——这是与标准 logit 的根本区别。

2.2 随机系数的分解

BLP 把随机系数进一步分解为"均值 + 人口统计 + 偏好冲击"三部分：

其中：

• $\beta$：所有消费者共同的系数均值（vector, dimension $K$）；
• $D_i$：消费者 $i$ 的可观察人口统计（income, age, family size…），dimension $D$；
• $\Pi$：$K \times D$ 矩阵，把人口统计映射到偏好；
• $v_i \sim N(0, I_K)$：消费者 $i$ 的不可观察偏好冲击（standard normal）；
• $\Sigma$：$K \times K$ 对角矩阵（实践中常取对角），控制偏好异质性的大小。

价格系数也可以写成 $\alpha_i = \alpha + \pi_y y_i + \sigma_\alpha v_{i,\alpha}$，其中 $y_i$ 是收入，体现"富人对价格不敏感"。

2.3 把效用拆成"共同部分 + 异质部分"

把 Equation 2.2 代入 2.1，效用可以拆成两块：

这个分解极其关键：

• $\delta_{jt}$ 只依赖产品 $j$ 的特征和均值系数，不随 $i$ 变化，叫 mean utility（平均效用）；
• $\mu_{ijt}$ 是消费者 $i$ 相对于均值的偏离，随 $i$ 变化，叫 consumer-product specific deviation；
• 外部产品（outside good）的效用归一化为 $u_{i0t} = \varepsilon_{i0t}$（即 $\delta_{0t}=0, \mu_{i0t}=0$）。

2.4 市场份额的表达式

由于 $\varepsilon_{ijt}$ 是 i.i.d. Extreme Value，条件在 $v_i, D_i$ 上，消费者 $i$ 选 $j$ 的概率仍是 logit 形式：

整个市场的份额 $s_{jt}$ 是个体选择概率在人口分布上的积分：

这里 $\theta_2 = \{\Pi, \Sigma, \pi_y, \sigma_\alpha\}$ 是"非线性参数"（nonlinear parameters）——它们出现在 $\mu_{ijt}$ 里，因为 $\mu$ 是 $v_i$ 的非线性函数。而 $\theta_1 = \{\beta, \alpha\}$ 是"线性参数"（linear parameters），只出现在 $\delta_{jt} = x_{jt}'\beta - \alpha p_{jt} + \xi_{jt}$ 中。

03 市场份额反演：收缩映射（Contraction Mapping）

BLP 估计的核心难点是：观测到的市场份额 $S_{jt}$（来自数据）和模型预测的份额 $s_{jt}(\delta, \theta_2)$（来自 Equation 2.5）之间有什么关系？我们要做的是"反演"——给定 $\theta_2$ 和观测份额 $S$，找到使 $s(\delta, \theta_2) = S$ 成立的 $\delta$。

3.1 为什么不能直接解 $\delta$？

在标准 logit 里，反演是解析的：

这是一个线性方程，直接 OLS/IV 即可。但在 BLP 中，Equation 2.5 是一个多维非线性积分，没有解析解。我们必须数值求 $\delta$。

3.2 BLP 收缩映射迭代

Berry (1994) 与 BLP (1995) 提出的迭代公式：

逐步解释：

• 给定初始猜测 $\delta^0$（通常设为 0，或用标准 logit 的 $\hat\delta$）；
• 用 Equation 2.5 在给定的 $\theta_2$ 下数值积分（用 Halton draws 或 pseudo-random draws 模拟 $F(D_i,v_i)$），算出模型预测份额 $s_{jt}(\delta^h, \theta_2)$；
• 比较观测份额 $S_{jt}$ 与预测份额 $s_{jt}$ 的对数差：如果 $S_{jt} > s_{jt}$，说明预测份额太低，需要提高产品 $j$ 的平均效用 $\delta_{jt}$，于是 $\delta^{h+1} > \delta^h$；
• 迭代直到 $\|\delta^{h+1} - \delta^h\|_\infty < \text{tol}$（例如 $10^{-8}$）。

3.3 为什么它收敛？（Blackwell 充分条件）

BLP 证明了这个映射是一个 contraction mapping：在标准条件下，存在 $\lambda \in (0,1)$ 使得 $\|g(\delta_1) - g(\delta_2)\| \leq \lambda \|\delta_1 - \delta_2\|$，其中 $g(\delta) = \delta + \ln S - \ln s(\delta)$。由 Banach 不动点定理，迭代收敛到唯一不动点 $\delta^*(\theta_2)$。

3.4 收敛后得到什么？

对给定的 $\theta_2$，收缩映射收敛后得到 $\delta^*(\theta_2)$。然后由 Equation 2.3 的定义：

注意：$\xi$ 既依赖 $\theta_2$（通过 $\delta^*$），也依赖 $\theta_1$（通过 $\beta, \alpha$）。但 $\theta_1$ 是线性的，我们可以在 GMM 内层解析地 concentrating out（见下一节）。

04 GMM 估计与工具变量

4.1 矩条件

$\xi_{jt}$ 是企业在定价时已知的未观察产品特征，因此 $p_{jt}$ 与 $\xi_{jt}$ 相关——这就是需求估计中的经典价格内生性（price endogeneity）。BLP 的识别依赖于工具变量：

其中 $Z_{jt}$ 是工具变量矩阵，dimension $L$（每行一个 moment）。

4.2 工具变量的构造

实践中常用三类 IV：

(a) Cost-side instruments（成本侧工具）

影响企业成本但不直接影响需求的变量：原材料价格、工资、汇率（进口车）。这些会通过边际成本影响价格，但与 $\xi_{jt}$（未观察需求冲击）无关。

(b) BLP instruments（基于产品特征的"竞争品特征和"）

对产品 $j$，构造：

直觉：如果产品 $j$ 所在市场有很多"类似特征"的竞争品，那么产品 $j$ 面临的竞争更激烈，均衡价格会更低；但这些竞争品的特征是外生的（由产品设计预先决定），与 $\xi_{jt}$ 不相关。在 multi-product firm 设定下，通常把"同企业其他产品特征"与"竞争企业产品特征"分开构造（Berry, Levinsohn, Pakes 1999 的最优 instruments）。

(c) Hausman instruments（其他市场的价格）

Hausman (1997) 提议用同一产品在其他市场的价格作为本市场价格的 IV：

逻辑：同一产品在其他市场的价格反映了共同的成本冲击（影响本市场价格），但与本市场的特定需求冲击 $\xi_{jt}$ 不相关。批评：如果全国性的品牌广告冲击同时影响所有市场，这个 IV 就不合法——实践中需要谨慎。

4.3 GMM 目标函数

定义堆叠后的矩向量：

GMM 目标函数：

其中 $W$ 是权重矩阵。两步 GMM：第一步用 $W = (Z'Z)^{-1}$ 得到 $\hat\theta^{(1)}$，第二步用 $W^* = (\frac{1}{NT}\sum Z_{jt}\xi_{jt}\xi_{jt}'Z_{jt}')^{-1}$（optimal GMM weight）得到 $\hat\theta^{GMM}$。

4.4 Concentrating out $\theta_1$（解析解）

由于 $\theta_1 = \{\beta, \alpha\}$ 只在 $\xi = \delta^*(\theta_2) - X\beta + \alpha p$ 中线性出现，对给定的 $\theta_2$，最优的 $\theta_1$ 可以解析地求出（类似 IV 回归的 closed form）：

因此外循环只需要在 $\theta_2$（通常不到 10 维）上优化，极大降低了维度灾难。

05 完整估计流程：内循环 + 外循环

BLP 估计是一个典型的嵌套优化（nested optimization）：

06 需求弹性与价格成本加成

估计完参数后，BLP 的两大应用是：(a) 计算 price elasticities；(b) 由 Bertrand-Nash 定价一阶条件反推出 marginal cost 与 markup。

6.1 Own-price elasticity

对个体 $i$，$P_{ijt}$ 对 $p_{jt}$ 的导数为 $\partial P_{ijt}/\partial p_{jt} = (-\alpha_i)(1 - P_{ijt}) P_{ijt}$。市场层面的弹性：

Cross-price elasticity：

关键性质：在 BLP 中，$\eta_{jk}$ 的大小取决于产品 $j$ 与 $k$ 在特征空间中的接近程度——相似产品之间 cross elasticity 大，这正是我们希望的。

6.2 Bertrand-Nash 一阶条件与 markup

假设企业 $f$ 同时拥有一组产品 $\mathcal{J}_f$，在差异化产品 Bertrand 竞争下，企业利润最大化：

一阶条件（matrix form）：

其中 $\Delta_{jk} = \partial s_j / \partial p_k$（$J \times J$ 矩阵），$\Omega_{jk} = 1$ 当 $j,k$ 同属一个企业 $f$，否则 0。解出 marginal cost：

这就是 BLP 在反垄断（antitrust）中被广泛使用的原因：估计出需求系统后，我们能反推出企业的真实边际成本与加价，进而模拟并购后的新均衡价格。

07 完整 Python 实现

下面用模拟数据演示 BLP 的核心算法：收缩映射 + GMM 外循环。为了可运行，我们简化为单一市场、只有价格与一个特征、随机系数只有价格的异质性 $\sigma_\alpha$。

"""
BLP (Berry-Levinsohn-Pakes) 简化实现
====================================
演示：单市场、J=10 个产品、R=100 个模拟消费者
关键组件：
  - 随机系数效用
  - 收缩映射（内循环）求 delta
  - GMM 外循环优化
依赖：numpy, scipy
"""
import numpy as np
from scipy.optimize import minimize

np.random.seed(42)

# ---------- 1. 模拟数据 ----------
J = 10           # 产品数
R = 100          # 模拟 draws 数
K = 2            # 特征维数（常数 + x1）

# 真实参数
beta_true = np.array([1.0, 0.5])     # 平均系数
alpha_true = -2.0                    # 平均价格系数
sigma_alpha_true = 1.0               # 价格系数的异质性

# 产品特征：常数项 + 一个观察特征
X = np.column_stack([np.ones(J), np.random.randn(J)])
p = np.exp(np.random.randn(J) * 0.3)  # 价格取正

# 未观察产品特征 xi
xi_true = 0.3 * np.random.randn(J)

# 模拟消费者的偏好冲击 v_i ~ N(0,1)，影响价格系数
v = np.random.randn(R)               # (R,)
# 消费者层面价格系数 alpha_i = alpha + sigma_alpha * v_i
alpha_i = alpha_true + sigma_alpha_true * v  # (R,)

# 真实 mean utility
delta_true = X @ beta_true + alpha_true * p + xi_true

# 计算真实市场份额（用于模拟"观测份额" S）
def calc_shares(delta, X, p, alpha_i):
    """
    输入：
      delta: (J,)  平均效用
      X:     (J,K) 特征
      p:     (J,)  价格
      alpha_i: (R,)  每个消费者的价格系数
    输出：
      s: (J,)  市场份额（在 R 个消费者上积分）
    """
    # 个体层面的 utility: (R, J)
    # mu_irj = - alpha_i_r * p_j  （这里只让价格系数异质）
    mu = -alpha_i[:, None] * p[None, :]      # (R, J)
    # 总 utility = delta_j + mu_irj
    u = delta[None, :] + mu                  # (R, J)
    # logit: P(r,j) = exp(u) / (1 + sum exp(u))
    u_max = u.max(axis=1, keepdims=True)
    exp_u = np.exp(u - u_max)                # (R, J)
    denom = 1.0 + exp_u.sum(axis=1, keepdims=True)  # (R,1)
    P = exp_u / denom                        # (R, J)
    # 在 R 个消费者上平均，得到市场份额
    s = P.mean(axis=0)                       # (J,)
    return s

S = calc_shares(delta_true, X, p, alpha_i)   # 观测份额（数据）

# ---------- 2. 工具变量 ----------
# 简化：用成本冲击作为 IV
cost_shock = 0.2 * np.random.randn(J)
# 构造 Z: 包含 X 和 cost_shock（与 p 相关、与 xi 独立）
Z = np.column_stack([X, cost_shock])          # (J, 3)

# ---------- 3. 收缩映射（内循环） ----------
def contraction_mapping(S_obs, X, p, sigma_alpha, v,
                       tol=1e-10, max_iter=1000):
    """
    给定 sigma_alpha，迭代求 delta 使得 s(delta) = S_obs
    返回收敛的 delta
    """
    delta = np.zeros(J)
    for h in range(max_iter):
        s_pred = calc_shares(delta, X, p, alpha_i=alpha_true + sigma_alpha * v)
        # 收缩映射：delta_new = delta + ln(S_obs) - ln(s_pred)
        delta_new = delta + np.log(S_obs) - np.log(s_pred)
        # 数值稳定：把 delta 减去均值（因为外部产品固定为 0）
        delta_new -= delta_new.mean()
        diff = np.max(np.abs(delta_new - delta))
        delta = delta_new
        if diff < tol:
            break
    return delta

# ---------- 4. GMM 目标函数（外循环） ----------
def gmm_objective(theta2, S_obs, X, p, Z, v):
    """
    theta2 = [sigma_alpha]  非线性参数
    其他参数 theta1 = [beta0, beta1, alpha] 通过 linear IV 解析集中化
    """
    sigma_alpha = theta2[0]
    # 内循环：求 delta
    delta = contraction_mapping(S_obs, X, p, sigma_alpha, v)
    # 集中化 theta1: 把 delta 对 [X, p] 用 IV Z 回归
    # Y = delta, W = [X, p], Z = instruments
    Wmat = np.column_stack([X, p])
    # IV 估计: theta1 = (W' Z (Z'Z)^-1 Z' W)^-1 W' Z (Z'Z)^-1 Z' delta
    ZtZ_inv = np.linalg.inv(Z.T @ Z)
    Pz = Z @ ZtZ_inv @ Z.T
    WtPzW = Wmat.T @ Pz @ Wmat
    WtPzd = Wmat.T @ Pz @ delta
    try:
        theta1 = np.linalg.solve(WtPzW, WtPzd)
    except np.linalg.LinAlgError:
        return 1e10
    # 结构误差 xi
    xi = delta - Wmat @ theta1
    # 矩条件
    g = Z.T @ xi / J                  # (L,)
    Wgmm = np.eye(Z.shape[1])         # 简化：第一步用单位权重
    obj = g @ Wgmm @ g
    return obj

# ---------- 5. 外循环优化 ----------
print("开始 GMM 优化 ...")
res = minimize(gmm_objective,
               x0=[0.5],
               args=(S, X, p, Z, v),
               method='L-BFGS-B',
               bounds=[(0.01, 5.0)],
               options={'maxiter': 50, 'ftol': 1e-8})
print(f"收敛状态: {res.message}")
print(f"估计 sigma_alpha = {res.x[0]:.3f}  (真实值 {sigma_alpha_true})")

# ---------- 6. 用估计的参数算弹性 ----------
sigma_alpha_hat = res.x[0]
delta_hat = contraction_mapping(S, X, p, sigma_alpha_hat, v)
alpha_i_hat = alpha_true + sigma_alpha_hat * v  # 注意：这里 alpha_true 简化了

# 计算 own-price elasticity
mu = -alpha_i_hat[:, None] * p[None, :]
u = delta_hat[None, :] + mu
u_max = u.max(axis=1, keepdims=True)
exp_u = np.exp(u - u_max)
denom = 1.0 + exp_u.sum(axis=1, keepdims=True)
P = exp_u / denom                  # (R, J)
s_hat = P.mean(axis=0)             # (J,)

# dP_j/dp_j = -alpha_i (1-P_j) P_j
dP_dp = -alpha_i_hat[:, None] * (1 - P) * P     # (R, J)
ds_dp = dP_dp.mean(axis=0)                       # (J,)
own_elasticity = (ds_dp * p / s_hat)
print("\nOwn-price elasticities (前5个产品):")
print(np.round(own_elasticity[:5], 3))

07+ 数据来源对照表（BLP 市场级数据）

BLP 用的不是个体级选择，而是市场级 aggregate 份额数据：每个市场 $t$ 中产品 $j$ 的份额 $S_{jt}$、价格 $p_{jt}$、特征 $x_{jt}$。下表对照经典与中国数据来源：

08 论文案例

09 常见错误与陷阱

10 进阶资料

• Conlon & Gortmaker (2020, JAE) — Best Practices for BLP Estimation. 必读，涵盖了过去 25 年 BLP 实现的所有坑和改进（加速、稳定、最优 IV）。
• Nevo (2000, JEL) — 上文已引，最实用的入门教程。
• Knittel & Metaxoglou (2014) — BLP 估计的数值敏感性：演示了不同优化器、tol、draws 选择如何改变点估计。
• PyBLP 开源包（jeffgortmaker/pyblp）— 现代 Python 实现，包含大量默认 best practices，可直接用于真实数据。
• Grigolon & Verboven (2014) — 带持久消费者需求（persistent consumer demand）的 BLP 扩展，处理耐用品市场。
• 毛亮、孙爽 等中文论文 — 用 BLP 估计中国汽车/乳制品/智能手机市场，是中文读者最好的对照参考。

方程总清单 / Equation Summary

本页全部方程按出现顺序汇总如下，共 23 个。每个方程均可在正文中找到对应的变量定义、设定理由与经济直觉。