教育与人力资本的结构估计：Mincer、Keane-Wolpin 与生命周期动态选择

"""
Keane-Wolpin 简化版 —— 3 期教育-工作选择模型反向归纳
=====================================================
状态:
  s ∈ {0,1}   教育水平 (0=高中, 1=大学)
  e ∈ {0,1,2} 工作经验
决策:
  d = 0 工作: 当期拿工资, s不变, e+1
  d = 1 上学: 当期付心理成本, s+1 (封顶1), e不变
工资方程 (Mincer):
  ln w = a0 + a1*s + a2*e + a3*e^2
终端价值 (t=3 之后):
  terminal(s,e) = gamma*s + delta*e   (退休/养老金)
冲击: 每期每选择 i.i.d. Type-I 极值 -> log-sum-exp
"""
import numpy as np
from scipy.special import logsumexp

# ---------- 1. 参数 ----------
BETA   = 0.9               # 贴现因子
A0     = 1.0               # 工资常数项
A1     = 0.15              # 教育回报率 (多一年/一级教育的工资增量)
A2     = 0.03              # 经验线性项
A3     = -0.002            # 经验二次项 (经验回报递减)
ETA    = 0.8               # 上学心理成本
GAMMA  = 3.0               # 终端价值中教育的权重
DELTA  = 1.0               # 终端价值中经验的权重
T      = 3                 # 总期数

# ---------- 2. 状态网格 ----------
S_GRID = [0, 1]            # 教育水平
E_GRID = [0, 1, 2]         # 经验
N_S, N_E = len(S_GRID), len(E_GRID)

def terminal_value(s, e):
    """退休后终端价值 (养老金/闲暇)"""
    return GAMMA * s + DELTA * e

def wage_util(s, e):
    """工作的当期效用 = Mincer 对数工资"""
    return A0 + A1 * s + A2 * e + A3 * e**2

# ---------- 3. 反向归纳 ----------
# V[t][s_idx, e_idx] = 期 t 开始时 (ex-ante) 的价值
# cv[t][d, s_idx, e_idx] = 选择 d 的 choice-specific value
V = {}
cv = {}

for t in range(T, 0, -1):
    V[t] = np.zeros((N_S, N_E))
    cv[t] = np.zeros((2, N_S, N_E))   # 2 个选择: 0=工作, 1=上学
    for si, s in enumerate(S_GRID):
        for ei, e in enumerate(E_GRID):
            # --- 选择 d=0: 工作 ---
            # 当期工资效用
            u_work = wage_util(s, e)
            # 下期状态: s 不变, e+1 (封顶)
            s_next, e_next = s, min(e + 1, 2)
            if t == T:
                fut_work = terminal_value(s_next, e_next)
            else:
                fut_work = V[t+1][S_GRID.index(s_next), E_GRID.index(e_next)]
            cv[t][0, si, ei] = u_work + BETA * fut_work

            # --- 选择 d=1: 上学 ---
            # 当期心理成本 (负效用)
            u_school = -ETA
            # 下期状态: s+1 (封顶), e 不变
            s_next, e_next = min(s + 1, 1), e
            if t == T:
                fut_school = terminal_value(s_next, e_next)
            else:
                fut_school = V[t+1][S_GRID.index(s_next), E_GRID.index(e_next)]
            cv[t][1, si, ei] = u_school + BETA * fut_school

            # ex-ante 价值 = logsumexp over d (EV1 冲击解析积分)
            V[t][si, ei] = logsumexp([cv[t][0, si, ei], cv[t][1, si, ei]])

# ---------- 4. 计算 CCP: 每期每状态选择"上学"的概率 ----------
def ccp_study(t, s, e):
    si, ei = s, e
    # P(study) = exp(v_school) / (exp(v_work)+exp(v_school))
    v_w = cv[t][0, si, ei]
    v_s = cv[t][1, si, ei]
    return np.exp(v_s) / (np.exp(v_w) + np.exp(v_s))

print("=== 各期各状态下'上学'的概率 P(study) ===")
print("      s=0,e=0  s=0,e=1  s=1,e=0  s=1,e=1")
for t in range(1, T+1):
    row = [ccp_study(t, s, e) for (s, e) in [(0,0),(0,1),(1,0),(1,1)]]
    print(f"t={t}:  " + "  ".join(f"{p:6.3f}" for p in row))

# ---------- 5. 模拟一个体的一生路径 ----------
def simulate_one(seed=42):
    rng = np.random.default_rng(seed)
    s, e = 0, 0     # 初始: 高中毕业, 无经验
    path = []
    for t in range(1, T+1):
        p_study = ccp_study(t, s, e)
        d = int(rng.random() < p_study)   # 抽选择
        path.append((t, s, e, d))
        # 状态转移
        if d == 0:   # 工作
            e = min(e + 1, 2)
        else:        # 上学
            s = min(s + 1, 1)
    return path

print("\n=== 模拟 20 个个体的教育-工作路径 ===")
print("(t, s, e, d=1表示上学)")
for i in range(5):
    path = simulate_one(seed=i)
    print(f"个体{i}: " + " -> ".join(f"(t{t},s{s},e{e},{'上学' if d==1 else '工作'})"
                                       for (t,s,e,d) in path))

# ---------- 6. 教育回报率的结构含义 ----------
# 比较: 若 A1=0 (教育无用), 与基准 A1=0.15 相比, 上学概率应大幅下降
print("\n=== 经济直觉 ===")
print(f"基准 A1={A1}: 第1期 s=0,e=0 时 P(study) = {ccp_study(1,0,0):.3f}")
print("若 A1 下降到 0.0 (教育无工资溢价), 上学动机应主要来自终端价值 GAMMA。")

*==============================================================
*  Mincer 工资方程: OLS / Heckman 选择修正 / 工具变量
*  数据: cfps_wage.dta (中国家庭追踪调查微观样本)
*  变量:
*    lnw        对数小时工资
*    edu        受教育年限
*    exp        潜在工作经验 (年龄-教育-6)
*    exp2       经验平方/100
*    participation  劳动参与哑变量 (1=工作)
*    college_proximity  Card(1995) 式IV: 成长地附近有无大学
*    mom_edu, father_income  家庭背景 (用于 Heckman 排除约束)
*==============================================================

clear all
set more off
use "cfps_wage.dta", clear

* --- (0) 构造变量 ---
gen exp2 = exp^2 / 100
lab var lnw            "对数小时工资"
lab var edu            "受教育年限"
lab var exp            "潜在经验"
lab var exp2           "经验平方/100"

* --- (1) OLS 基准 Mincer 方程 (含能力偏差, 仅作对照) ---
reg lnw edu exp exp2 i.province i.year, robust
estimates store OLS
* 系数 edu 是"简单教育回报率", 但混杂能力偏差

* --- (2) Heckman 选择模型: 修正"只有工作者才观察到工资"的选择偏误 ---
* 第一阶段: 劳动参与方程 (排除变量 mom_edu, father_income 不进工资方程)
* 第二阶段: 工资方程 (加逆 Mills ratio lambda)
heckman lnw exp exp2 i.province i.year, ///
    select(participation = edu mom_edu father_income exp exp2) ///
    robust
estimates store Heckman
* 若 lambda 显著, 说明存在选择偏误, Heckman 修正必要

* --- (3) 工具变量 IV: Card(1995) college proximity 作教育的 IV ---
* 排他性约束: 成长地附近有无大学影响上学成本 (相关性),
*             但不直接影响成年后的工资 (外生性)
ivregress 2sls lnw exp exp2 i.province i.year ///
    (edu = college_proximity), robust
estimates store IV
estat firststage          // 检查弱工具变量 (F>10)
estat overid              // 过度识别检验 (若有多个IV)

* --- (4) 对比三种估计的教育回报率 ---
estimates table OLS Heckman IV, ///
    b(%9.4f) star(0.1 0.05 0.01) keep(edu)
* 典型结果: OLS 偏低(能力偏差方向取决于假设),
*           Heckman/IV 通常给出更高的"边际"教育回报率