MPEC 方法：带均衡约束的数学规划

01 概念直觉：MPEC 到底改了什么

传统结构估计（如 Rust 的 NFXP, Nested Fixed Point）是一种双层嵌套算法：外层优化器不断猜测参数 $\theta$；对每一个 $\theta$，内层都要先"求解模型均衡"（迭代 Bellman 方程求价值函数 $V(\theta)$，或不动点求市场份额 $s(\theta)$），再把解得的 $V$ 塞进似然/矩目标函数，比较是否更优。这种"每试一个参数就重解一次模型"的做法，计算成本随维度爆炸——外层次数 × 内层不动点迭代次数 × 每次迭代的模拟量。

MPEC 的核心思想：干脆不要内循环了。把价值函数 $V$（或市场份额 $\delta$）当作与 $\theta$ 并列的决策变量，把"$V$ 必须满足 Bellman 均衡条件"写成一个等式约束 $g(\theta,V)=0$。于是整个估计问题变成一个标准的、单层的、带约束的非线性规划（NLP）：

经济学直觉：我们真正要找的是"一组同时满足两个条件的解"——① 它让似然/矩目标最大（拟合数据）；② 它满足模型的均衡方程（个体最优化、市场出清）。MPEC 只是老实地把这两个条件同时写进优化问题，让现代化的内点法 NLP 求解器去处理它。理论上，这与 NFXP 估计的是同一个解；但计算上，它消除了"每步重解不动点"的浪费。

02 MPEC 的一般形式

设 $\theta$ 为结构参数向量，$V$ 为模型内生均衡变量向量（在动态离散选择中是价值函数，在 BLP 中是产品的 mean-utility $\delta_j$），$f(\theta,V)$ 为目标函数（负对数似然或 GMM 二次型）。MPEC 的标准形式为：

其中 $g(\theta,V)=0$ 就是均衡约束——它强制"价值函数必须满足 Bellman 方程"或"市场份额必须等于需求预测份额"。求解时，$\theta$ 和 $V$ 是同步被调整的：求解器在降低 $f$ 的同时，必须保证约束残差 $g$ 趋近于 0。由于不强制每步都精确解出不动点，内点法可以沿着"约束越来越紧、目标越来越好"的路径大步前进，这正是它比嵌套法快的原因。

03 MPEC vs NFXP：计算效率与收敛性对比

04 动态离散选择：Rust 公车问题的 MPEC 形式

Rust (1987) 的公车更换模型中，状态为里程 $x$，司机选择"更换"($i=1$)或"继续运营"($i=0$)。Bellman 方程为

在 NFXP 里，每猜一次 $\theta=(RC,c_0,c_1,\beta)$，都要把上式迭代到收敛得到 $V(\theta)$。MPEC 改写：把每个状态点的 $V(x)$ 都当作变量，Bellman 方程从"必须恒等式"变成一个约束——注意 logit 极值误差下，Bellman 方程有解析的平滑形式（用 logsum 替换 max）：

目标函数是数据的（条件选择概率）负对数似然。于是整个问题就是：在"每个状态的 $V(x)$ 必须满足平滑 Bellman 残差为 0"的约束下，最小化负对数似然。变量是 $(\theta, \{V(x)\}_{x=1}^N)$，约束数等于状态点数 $N$。这正是 Su & Judd (2012) 对 Rust 问题的 MPEC 重构，报告显示比 NFXP 快得多且更稳定。

05 约束条件的构造：Bellman 约束 / 市场份额约束

5.1 Bellman 方程约束（动态离散选择）

约束就是"Bellman 残差"。对每个状态 $x$、每个决策时点，写出 $g_x(\theta,V)=V(x)-\text{logsum}(\theta,V)=0$。构造要点：

• 用 logsum 平滑版而非硬 max，否则约束不可微，NLP 求解器无法求 Jacobian。
• 约束数 = 状态空间维数 $N$，与 $\theta$ 分离，稀疏性强——这正是 KNITRO/IPOPT 擅长的结构。

5.2 市场份额约束（BLP / 需求估计）

在 BLP 中，BLP 不动点要解出 mean-utility $\delta$ 使预测份额等于观测份额。MPEC 把 $\delta$ 当变量，市场出清直接写成约束：

即"模型预测的第 $j$ 个产品份额"减去"观测份额"必须为 0。$s_j^{model}$ 由随机系数 logit 需求积分给出。这样就完全省去了 BLP 内嵌的不动点循环（经典 BLP 每估计一步要解 $\delta$ 的不动点，极易发散）。

06 BLP 中的 MPEC：Dubé-Fox-Su (2012)

BLP (1995) 的经典估计是"三层嵌套"：① 解 $\delta$ 不动点；② 用 BLP 工具变量做价格内生性的矩条件；③ 外层 GMM 优化随机系数参数。Dubé, Fox & Su (2012, RAND Journal) 指出这种嵌套在大规模市场下计算上极不稳定：内层不动点的误差会污染外层 GMM 梯度，导致估计发散或收敛到错误值。

他们把整个 BLP 重写为一个 MPEC：变量 = (随机系数参数 $\theta^p$, 产品 mean-utility $\delta_{jt}$)；约束 = 市场份额出清 $s_{jt}^{model}(\delta,\theta^p)=s_{jt}^{obs}$；目标 = BLP GMM 二次型（含工具变量矩条件）。并用商用内点求解器 KNITRO、开源 IPOPT 直接求解。结论是：MPEC 化的 BLP 更稳、更快、数值上更可信，成为后续大规模需求估计的标准做法。

07 完整 Python 代码：scipy 约束优化实现 MPEC

下面用一个可直接运行的最小 MPEC 实例：一个 3 状态的 logit 动态离散选择模型。把 $\theta=(RC,c_1)$ 与三个价值函数 $V=(V_1,V_2,V_3)$ 一起作为变量，Bellman 残差作为等式约束，最小化一个简单的"矩/似然"目标。用 scipy.optimize.minimize(method='SLSQP') 实现。

import numpy as np
from scipy.optimize import minimize

# ============================================================
# MPEC 完整实例：3状态 logit 动态离散选择（Rust 风格）
# 决策：i=0 继续, i=1 更换。固定 beta=0.95
# 状态 x=1,2,3，运营成本 c(x)=c1*x；更换后回到 x=1
# 决策变量 z = [theta=(RC, c1), V=(V1,V2,V3)]，共 5 维
# 约束：3 个 Bellman 残差 g_k = V_k - logsum(theta,V) = 0
# 目标：让"模型选择概率"贴近一组给定的"观测"选择概率
# ============================================================

beta = 0.95
# 观测到的"继续"概率（模拟一份数据矩/选择频率，用来拟合理想）
obs_p_continue = np.array([0.90, 0.60, 0.30])   # x越大越倾向更换

def logsum(a, b):
    """log(exp(a)+exp(b)) 的数值稳定写法（Bellman 平滑算子）"""
    m = max(a, b)
    return m + np.log(np.exp(a - m) + np.exp(b - m))

def unpack(z):
    RC, c1 = z[0], z[1]          # 结构参数
    V = np.array([z[2], z[3], z[4]])  # 价值函数（与 theta 并列的变量）
    return RC, c1, V

def bellman_constraints(z):
    """MPEC 的核心：3 个 Bellman 残差约束 g_k(theta,V)=0"""
    RC, c1, V = unpack(z)
    g = np.zeros(3)
    for k, x in enumerate([1, 2, 3]):
        # 选项1：更换 -> 回到状态1，成本 RC
        val_replace = -RC + beta * V[0]
        # 选项0：继续 -> 状态 x+1（x=3 时留在3），成本 c1*x
        x_next = min(x + 1, 3)
        val_operate = -c1 * x + beta * V[x_next - 1]
        # Bellman 残差 = V_k - logsum(两选项)
        g[k] = V[k] - logsum(val_replace, val_operate)
    return g

def choice_probs(z):
    """在给定 (theta,V) 下，导出各状态的'继续'概率（logit）"""
    RC, c1, V = unpack(z)
    p = np.zeros(3)
    for k, x in enumerate([1, 2, 3]):
        val_replace = -RC + beta * V[0]
        x_next = min(x + 1, 3)
        val_operate = -c1 * x + beta * V[x_next - 1]
        # 继续概率 = softmax(operate, replace)
        p[k] = np.exp(val_operate) / (np.exp(val_operate) + np.exp(val_replace))
    return p

def objective(z):
    """目标：模型继续概率与观测概率的加权距离（GMM/似然式）"""
    p = choice_probs(z)
    return np.sum((p - obs_p_continue)**2)   # 二次型矩距离

# ---------- 组装 MPEC 问题 ----------
# 变量初值：theta 猜 (RC=5, c1=1)，V 猜 0
z0 = np.array([5.0, 1.0, 0.0, 0.0, 0.0])

# 等式约束：Bellman 残差 = 0 （3 个）
cons = ({'type': 'eq', 'fun': bellman_constraints})

# 参数与 V 的合理边界（RC>0, c1>0）
bounds = [(0.01, 50), (0.01, 20), (-1e3, 1e3), (-1e3, 1e3), (-1e3, 1e3)]

# 用 SLSQP 求解（小规模问题足够；大规模换 cyipopt/KNITRO）
res = minimize(objective, z0, method='SLSQP', bounds=bounds,
               constraints=cons,
               options={'maxiter': 500, 'ftol': 1e-10, 'disp': True})

RC_hat, c1_hat, V_hat = unpack(res.x)
print("\n===== MPEC 估计结果 =====")
print("收敛状态:", res.message)
print("RC_hat = %.4f, c1_hat = %.4f" % (RC_hat, c1_hat))
print("V_hat =", np.round(V_hat, 4))
print("Bellman 残差（应≈0）:", np.round(bellman_constraints(res.x), 6))
print("模型继续概率:", np.round(choice_probs(res.x), 4),
      " 观测:", obs_p_continue)

08 论文案例

09 MPEC 常见陷阱（≥3）与进阶资料

进阶资料

• Su & Judd (2012) 原文附录：MPEC 在 Rust 问题上的完整约束构造。
• Dubé, Fox, Su (2012)：BLP-MPEC 的 KNITRO/IPOPT 实现细节与诊断。
• Nocedal & Wright, Numerical Optimization Ch.19（内点法/约束优化理论）。
• pyblp / KyPlus 等开源 BLP 包：BLP 的 MPEC 化工程参考。