生产函数估计：Olley-Pakes、Levinsohn-Petrin 与 Ackerberg-Caves-Frazer

01 经济环境与假设

生产函数估计是产业组织、宏观与贸易实证中最基础也最棘手的问题。我们关心企业层面的全要素生产率（TFP）以及资本、劳动、中间投入的产出弹性。考虑一个 Cobb-Douglas 生产函数（对数形式）：

其中 $y_{it}$ 是企业 $i$ 在 $t$ 期的产出（对数），$k_{it}, l_{it}, m_{it}$ 分别是资本、劳动、中间投入（均取对数）。误差项被分解为两部分：$\omega_{it}$ 是企业可观测但计量经济学家不可观测的生产率冲击，$\eta_{it}$ 是独立同分布的测量误差或意外冲击。我们的目标是一致估计 $\beta_k, \beta_l, \beta_m$，再用 $\hat{\omega}_{it}=y_{it}-\hat{\beta}_0-\hat{\beta}_k k_{it}-\hat{\beta}_l l_{it}-\hat{\beta}_m m_{it}$ 作为 TFP。

1.1 市场结构假设：完全竞争 vs 垄断竞争

OP/LP/ACF 框架默认企业在产品市场上是价格接受者或面对给定需求曲线，产出 $Y_{it}$ 用实物量或经行业价格平减的销售额衡量。若企业有市场势力（markup $\mu > 1$），则用收入产出 $r_{it}=p_{it}+y_{it}$ 估计会把 markup 混进 TFP，得到的是 revenue-TFP 而非 physical-TFP。Foster-Haltiwanger-Syverson (2008) 区分了这两者：在同质性行业中两者接近，在差异化行业中 revenue-TFP 与 markup 共变，后续回归需谨慎。劳动市场同理：默认工资外生给定；若企业有买方势力，$\beta_l$ 的估计也会失真。

1.2 要素调整成本假设：自由要素与准固定要素

这是整个框架的关键时序假设。劳动 $l_{it}$ 与中间投入 $m_{it}$ 是"自由"或"可变"要素：企业在 $t$ 期观察到 $\omega_{it}$ 后即可无成本调整，因此 $l_{it}, m_{it}$ 与 $\omega_{it}$ 同期相关。资本 $k_{it}$ 是"准固定"要素：它由前期投资累积，$k_{it}=(1-\delta)k_{i,t-1}+i_{i,t-1}$，在 $t$ 期期初就已确定，企业在选择 $k_{it}$ 时只能看到 $t-1$ 期及以前的生产率信息。这一"资本前定"假设是后文矩条件 $E[\xi_{it}\,|\,k_{it}]=0$ 的根基——资本对当期生产率创新正交，但对当期 $\omega_{it}$ 本身不正交（因为资本累积基于预期）。

1.3 生产率过程假设：一阶马尔可夫

假设生产率遵循一阶马尔可夫过程：

$\xi_{it}$ 是当期"未预期到"的生产率创新。企业在 $t-1$ 期做投资决策时只知道 $\omega_{i,t-1}$，不知道 $\xi_{it}$，所以 $k_{it}$ 与 $\xi_{it}$ 正交。这给了我们识别 $\beta_k$ 的外生变化。若马尔可夫阶数错了（如实际是 AR(2)），矩条件会失效。

1.4 企业进入与退出假设

企业在每期期初决定是否继续经营。设保留值 $\underline{\omega}_{it}(k_{it})$，当 $\omega_{it} \ge \underline{\omega}_{it}(k_{it})$ 时存活（$\chi_{it}=1$），否则退出。资本越多的企业"留下"所需的最低生产率越低（资本是沉没缓冲），因此在存活样本中 $\omega_{it}$ 的分布是左截断的，且截断点随 $k_{it}$ 变化。这就是选择偏误的来源。

02 内生性：同时性偏误与选择偏误

2.1 同时性偏差 (Simultaneity Bias)

企业在每期知道自己的 $\omega_{it}$，并据此选择当期劳动 $l_{it}$ 与中间投入 $m_{it}$。因此 $\omega_{it}$ 与 $l_{it}, m_{it}$ 相关：

用 OLS 直接回归，劳动与中间投入系数被向上偏误：高生产率企业多雇劳动，OLS 把"生产率带来的高产出"错记为"劳动的高边际产出"。资本 $k_{it}$ 因前定而内生性较弱。

2.2 选择偏差 (Selection / Attrition Bias)

低生产率企业更可能退出。若样本只含存活企业，则给定资本下 $\omega_{it}$ 分布左截断：

这种截断使 $k_{it}$ 与残差负相关，OLS 低估资本系数 $\beta_k$。OP 用美国电信设备业数据发现：OLS 估计的规模报酬远小于 1，修正后接近常数——差别主要来自选择偏误。

03 Olley-Pakes (1996)：投资作为代理变量

3.1 单调性假设

OP 的关键假设：企业最优投资 $i_{it}$ 是当期 $\omega_{it}$ 与资本 $k_{it}$ 的严格递增函数（控制 $k$ 后）：

严格单调性要求 $i_{it}>0$（投资为零时逆函数不存在）。代入生产函数得：

3.2 三步法

04 Levinsohn-Petrin (2003)：中间投入代理

LP 用中间投入 $m_{it}$ 代替投资作为代理变量：

经济学理由：中间投入几乎恒为正，避免 OP 零投资丢样本；中间投入调整成本低、对 $\omega$ 反应更连续平滑。三步法结构与 OP 平行：

• 第一阶段：用 $(m_{it},k_{it})$ 的四阶多项式非参逼近 $\phi(m_{it},k_{it})$，OLS 得 $\hat{\beta}_l$。
• 第二阶段：对候选 $(\beta_k^*,\beta_m^*)$，构造 $\hat{\omega}_{it}=\hat{\phi}_{it}-\beta_k^* k_{it}-\beta_m^* m_{it}$，对 $\omega_{i,t-1}$ 非参回归得 $E[\omega_{it}\mid\omega_{i,t-1}]$，从而得创新 $\xi_{it}(\beta^*)$。
• 第三阶段：用矩 $E[\xi_{it}\cdot k_{it}]=0$ 与 $E[\xi_{it}\cdot m_{i,t-1}]=0$ 最小化，得 $\hat{\beta}_k,\hat{\beta}_m$。

05 Ackerberg-Caves-Frazer (2015)：函数相关性批判

ACF 指出 OP/LP 第一阶段存在根本的函数相关性 (functional dependence)：若 $l_{it}$ 也是 $(\omega_{it},k_{it})$ 的函数 $l_{it}=l(\omega_{it},k_{it})$，则 $\beta_l l_{it}$ 会被 $\phi(m_{it},k_{it})$ 完全吸收，第一阶段 $\beta_l$ 根本不可识别。

5.1 Value-added vs Gross Output

ACF 强调：若生产函数是 gross output（总产出），必须把 $m_{it}$ 放进生产函数本身，不能先减 $m$ 得 value-added 再估计——因为 $m_{it}$ 既是要素又是代理变量，两步处理会出错。推荐直接估计 gross output 生产函数。

06 识别条件与矩条件

整套方法的识别依赖三个互相关联的假设，缺一不可：

1. 单调性：$i(\omega,k)$ 或 $m(\omega,k)$ 在控制 $k$ 后关于 $\omega$ 严格单调，逆函数存在。若企业存在产能约束、投资不可逆，单调性会被破坏。
2. 当期自由/前定要素的时序：$l_{it},m_{it}$ 同期响应 $\omega_{it}$；$k_{it}$ 前定于 $\omega_{it}$ 的创新 $\xi_{it}$。这是 $E[\xi_{it}\cdot k_{it}]=0$ 的来源。
3. 一阶马尔可夫：$\omega_{it}$ 的条件期望只依赖 $\omega_{i,t-1}$，滞后二阶及以上与创新正交，提供了 $l_{i,t-1},m_{i,t-1}$ 作为额外工具。

6.1 生产率创新的构造（补全）

6.2 各方法矩条件（补全，禁止"常用矩汇总"跳步）

常用矩条件汇总：$E[\xi_{it}\mid k_{it}]=0$、$E[\xi_{it}\mid l_{i,t-1}]=0$、$E[\xi_{it}\mid m_{i,t-1}]=0$、$E[\xi_{it}\mid k_{i,t-1}]=0$。过度识别检验（Hansen J）可用来检验工具外生性。Wooldridge (2009) 把这些矩写成单步 GMM，标准误由 sandwich 矩阵直接给出，避免两步渐近修正。

07 数据要求与来源

估计 OP/LP/ACF 需要企业面板数据，至少包含：

常用数据源：中国工业企业数据库（1998–2007，规模以上工业企业，杨汝岱、聂辉华等大量研究基于此）、美国 Census LRD/ASM（OP、LP 原始论文所用）、智利/哥伦比亚制造业面板（LP 原文用智利）。中国数据需注意：2008 年后工业企业库口径变化、规模以下企业"消失"不等于退出，退出定义要谨慎；资本存量用永续盘存法 $K_{it}=(1-\delta)K_{i,t-1}+I_{it}$ 构造。此外，OP 要求 $i_{it}>0$ 的样本截断会损失约 15%–30% 的观测，LP 用中间投入则几乎不损失样本——这一点在比较方法时必须同时报告样本量。

08 应用解读：TFP 与后续分析

估计出 $\hat{\beta}$ 后，TFP 残差 $\hat{\omega}_{it}$ 是后续研究的核心变量。常见用法：

• TFP 增长分解：Olley-Pakes 分解 $\Delta \bar{\omega}=\Delta \bar{\omega}_{unweighted}+\Delta \sum(s_{it}-\bar{s})(\omega_{it}-\bar{\omega})$，把总量 TFP 增长拆成"企业自身进步"与"资源再配置 (reallocation)"。中国制造业研究普遍发现 reallocation 贡献远大于平均进步。
• 与出口/创新/外资的关系：用 $\hat{\omega}_{it}$ 对出口虚拟变量、专利、外资份额回归，考察"自我选择"（高效率企业出口）与"学习效应"（出口提升效率）。
• 资源误置：Hsieh-Klenow (2009) 框架用 $\hat{\omega}$ 的离散度衡量 misallocation，反事实"如果消除误置 TFP 能提升多少"。
• 政策评估：用 $\hat{\omega}_{it}$ 作为企业绩效变量，评估国企改革、开发区、减税政策的因果效应。
• 跨国/跨行业比较：把 $\hat{\omega}$ 标准化（行业-年份均值归零）后比较中国与美国、印度的 TFP 分布离散度，是 Hsieh-Klenow 框架的标准做法。

09 完整代码：Stata + Python

9.1 Stata：opreg / levpet / acfest

*==============================================================*
* 生产函数估计：OP / LP / ACF 完整示例
* 数据：企业面板 id-year，变量 y(增加值/总产出) k l m
* 安装：ssc install opreg, replace
*       ssc install levpet, replace
*       ssc install acfest, replace
*==============================================================*
clear all
set more off
use "firm_panel.dta", clear

xtset id year
gen ln_y  = ln(value_added)
gen ln_k  = ln(capital)
gen ln_l  = ln(labor)
gen ln_m  = ln(material)
gen ln_inv= ln(investment)

* --- 0. OLS 参照 ---
reg ln_y ln_k ln_l, robust
estimates store ols

* --- 1. Olley-Pakes (1996) ---
opreg ln_y = ln_l, state(ln_k) proxy(ln_inv) exit(exit_flag) ///
    vce(boot, reps(50))
estimates store op

* --- 2. Levinsohn-Petrin (2003) ---
levpet ln_y, free(ln_l) state(ln_k) proxy(ln_m) valueadded
estimates store lp_va

* gross output 形式
levpet ln_y_gross, free(ln_l) state(ln_k) proxy(ln_m)
estimates store lp_go

* --- 3. Ackerberg-Caves-Frazer (2015) ---
acfest ln_y_gross, free(ln_l) state(ln_k) proxy(ln_m) ///
    vce(boot, reps(100)) tech(nr)
estimates store acf

* --- 4. 导出 TFP 对比 ---
predict tfp_ols, xb
gen tfp_op = ln_y - _b[ln_l]*ln_l - _b[ln_k]*ln_k
predict tfp_acf, omega
esttab ols op lp_va acf, se star(* 0.1 ** 0.05 *** 0.01) ///
    mtitles("OLS" "OP" "LP(VA)" "ACF")

9.2 Python：ACF 三步法手写实现

# ============================================================
# ACF (2015) 生产函数估计的最小可运行实现
# 第一阶段: (m,k) 多项式拟合 phi_hat (不含 l)
# 第二/三阶段: GMM 同时估计 (beta_l, beta_k, beta_m)
# ============================================================
import numpy as np
import pandas as pd
from scipy.optimize import minimize
import statsmodels.api as sm

def acf_estimate(df, y='y', l='l', k='k', m='m', id_col='id', t_col='t'):
    df = df.sort_values([id_col, t_col]).reset_index(drop=True)
    # 第一阶段: phi(m,k) = beta0 + h*(m,k)，排除 l
    Xp = _poly2d(df[m].values, df[k].values, degree=4)
    Xp = sm.add_constant(Xp)
    df['phi_hat'] = sm.OLS(df[y].values, Xp).fit().fittedvalues

    def obj(params, df):
        bl, bk, bm = params
        df = df.copy()
        df['w'] = df['phi_hat'] - bl*df[l] - bk*df[k] - bm*df[m]
        df['wl'] = df.groupby(id_col)['w'].shift(1)
        v = df.dropna(subset=['wl']).copy()
        x = v['wl'].values; yv = v['w'].values
        Xq = np.column_stack([np.ones_like(x), x, x**2, x**3])
        c, *_ = np.linalg.lstsq(Xq, yv, rcond=None)
        xi = yv - Xq @ c
        kt = v[k].values
        ll = v.groupby(id_col)[l].shift(1).reindex(v.index).values
        ml = v.groupby(id_col)[m].shift(1).reindex(v.index).values
        Z = np.column_stack([np.ones_like(kt), kt, ll, ml])
        g = np.nanmean(xi[:, None]*Z, axis=0)
        return g @ np.eye(len(g)) @ g

    init = sm.OLS(df[y], sm.add_constant(df[[l,k,m]])).fit().params
    x0 = [init[l], init[k], init[m]]
    r = minimize(obj, x0, args=(df,), method='L-BFGS-B',
                 bounds=[(0,1)]*3)
    bl, bk, bm = r.x
    df['tfp'] = df[y] - bl*df[l] - bk*df[k] - bm*df[m]
    return dict(beta_l=bl, beta_k=bk, beta_m=bm, tfp=df['tfp'].values)

def _poly2d(x1, x2, degree=4):
    terms = []
    for i in range(1, degree+1):
        for j in range(i+1):
            terms.append((x1**j)*(x2**(i-j)))
    return np.column_stack(terms)

# ---- 模拟数据: 1000 企业 x 10 年 ----
np.random.seed(42); N, T = 1000, 10
ids = np.repeat(np.arange(N), T); ts = np.tile(np.arange(1,T+1), N)
bl, bk, bm = 0.60, 0.20, 0.15
omega = np.zeros(N*T)
for t in range(1, T):
    omega[t::T] = 0.8*omega[t-1::T] + np.random.normal(0,0.1,N)
k = np.log1p(np.cumsum(np.random.exponential(0.1,N*T),axis=0)/T)
l = 0.5*omega + np.random.normal(0,0.1,N*T)
m = 0.3*omega + 0.2*k + np.random.normal(0,0.05,N*T)
y = bl*l + bk*k + bm*m + omega + np.random.normal(0,0.02,N*T)
d = pd.DataFrame({'id':ids,'t':ts,'y':y,'l':l,'k':k,'m':m})
res = acf_estimate(d)
print(f"bl={res['beta_l']:.3f}(0.60) bk={res['beta_k']:.3f}(0.20) bm={res['beta_m']:.3f}(0.15)")

10 论文案例与常见错误

10.1 经典论文

10.2 常见错误

方程总清单 / Equation Summary

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