生命周期消费-储蓄结构模型：欧拉方程、借贷约束与缓冲存货

# -*- coding: utf-8 -*-
# 代码块1: EGM 求解生命周期消费函数 (Carroll 2006 / Deaton 1991)
import numpy as np
from numpy.polynomial.hermite_e import hermegauss

# ---------- 1. 模型参数 ----------
T, R   = 50, 40        # 总期数 / 工作期数 (后10期退休)
beta   = 0.96          # 主观贴现率 (impatient: beta(1+r)<1)
r      = 0.02          # 无风险利率
Rf     = 1.0 + r
sigma  = 2.0           # 相对风险厌恶 (EIS=1/sigma)
sig_eta, sig_eps = 0.10, 0.10   # 持久/暂时冲击标准差 (log)
a_min, y_R = 0.0, 0.60          # 借贷约束 / 退休养老金
N_a, N_quad = 60, 7             # 资产网格数 / 求积点数

# 收入冲击 Gauss-Hermite 求积 (对数正态)
nodes, w = hermegauss(N_quad); w = w/np.sqrt(2*np.pi)
eta_grid = np.exp(sig_eta*np.sqrt(2)*nodes)   # 持久冲击乘子
eps_grid = np.exp(sig_eps*np.sqrt(2)*nodes)   # 暂时冲击乘子
up    = lambda c: c**(-sigma)                 # 边际效用
invup = lambda v: v**(-1.0/sigma)             # 边际效用逆函数

# ---------- 2. 标准化: 以持久收入 p_t 为单位 ----------
# m_t = x_t/p_t 归一化现金, C_t(m) 归一化消费函数
a_next_grid = np.linspace(a_min, 20.0, N_a)   # 外生下一期资产网格
policy = {T: lambda m: np.clip(m, 1e-10, None)}   # 终端: 吃光 c_T=x_T

# ---------- 3a. 退休期反向归纳 (无冲击, 解析年金解) ----------
for t in range(T-1, R, -1):
    k = T - t + 1                              # 含本期的剩余期数
    disc = (beta*Rf)**(1.0/sigma)
    ann  = sum([(disc/Rf)**j for j in range(k)])  # 年金现值因子
    policy[t] = lambda m, ann=ann: np.clip(m/ann, 1e-10, m)

# ---------- 3b. 工作期 EGM (有收入冲击) ----------
for t in range(R, 0, -1):
    pn = policy[t+1]                           # 下一期消费函数
    me = np.zeros(N_a); ce = np.zeros(N_a)
    for i, ap in enumerate(a_next_grid):        # 固定下一期资产 a'
        rhs = 0.0
        for ie, eta in enumerate(eta_grid):    # 持久冲击
            for je, eps in enumerate(eps_grid): # 暂时冲击
                mn = Rf*ap/eta + eps            # 下一期归一化现金
                cn = pn(mn)                     # 下一期归一化消费
                rhs += w[ie]*w[je]*up(eta*cn)   # E_t[u'(c_{t+1})]
        c_t = invup(beta*Rf*rhs)                # 解析解本期消费
        m_t = c_t + ap                          # 内生现金网格
        me[i], ce[i] = m_t, c_t
    # 借贷约束: 最小现金处 c=m (即 a=a_min=0)
    if me[0] > a_min:
        me = np.concatenate([[a_min], me]); ce = np.concatenate([[a_min], ce])
    o = np.argsort(me); me, ce = me[o], ce[o]  # 单调化
    policy[t] = lambda m, me=me, ce=ce: np.interp(m, me, ce, left=me[0], right=me[-1])

# ---------- 4. 输出检验 ----------
print("工作期 t=1 的消费函数 C(m):")
for m in [1.0, 3.0, 6.0, 10.0]:
    print(f"  m={m:5.1f}  C={policy[1](m):.4f}  储蓄率={1-policy[1](m)/m:.3f}")

# -*- coding: utf-8 -*-
# 代码块2: 模拟生命周期路径 + 缓冲存货目标资产
# (接续代码块1的 policy / 参数; 此处给出完整模拟部分)
import numpy as np
rng = np.random.default_rng(2024)
N = 5000                                   # 模拟家庭数
p = np.ones(N); a = np.zeros(N)           # 初始: 持久收入=1, 资产=0
C = np.zeros((N,T)); A = np.zeros((N,T)); Y = np.zeros((N,T))
for t in range(1, T+1):
    if t <= R:                             # 工作期
        eta = np.exp(rng.normal(0, sig_eta, N))
        eps = np.exp(rng.normal(0, sig_eps, N))
        p = p*eta; y = p*eps
    else:                                  # 退休期: 固定养老金
        y = np.full(N, y_R*np.median(p))
    x = Rf*a + y                           # 现金在手
    c = p*policy[t](np.clip(x/p, 0, 20))   # 最优消费
    a = np.maximum(x - c, a_min)           # 期末资产(不违反约束)
    C[:,t-1], A[:,t-1], Y[:,t-1] = c, a, y

# 跨家庭平均生命周期曲线
print("年龄  平均收入  平均消费  平均资产")
for t in [0,9,19,29,39,44,49]:
    print(f" {t+1:2d}   {Y[:,t].mean():6.3f}  {C[:,t].mean():6.3f}  {A[:,t].mean():7.3f}")

# 缓冲存货目标: 年轻工作期资产/收入比的横截面均值
ratio = np.nanmean(A[:,4:25]/np.maximum(Y[:,4:25], 1e-6))
print(f"\n年轻工作期(年龄5-25)平均资产/收入比 ≈ {ratio:.3f}")

# -*- coding: utf-8 -*-
# 代码块3: 欧拉方程 GMM 估计 (线性化, 滞后变量作工具)
import numpy as np
from scipy.optimize import minimize
lnC = np.log(np.maximum(C, 1e-6))
dlnC = lnC[:,1:] - lnC[:,:-1]             # Delta ln c_{t+1}
lnY = np.log(np.maximum(Y, 1e-6))
dlnY = lnY[:,1:] - lnY[:,:-1]

# 模型: Delta ln c_{t+1} = delta0 + eps,  工具 z = [滞后消费增长, 滞后收入增长]
dep = dlnC[:,1:].ravel()                   # 被解释变量
z1  = dlnC[:,:-1].ravel()                  # 滞后一期消费增长
z2  = dlnY[:,:-1].ravel()                  # 滞后一期收入增长
Z   = np.column_stack([np.ones_like(dep), z1, z2])

def gmm_obj(theta):
    e = dep - theta[0]                     # 欧拉残差
    g = (Z * e[:,None]).mean(axis=0)       # 矩条件
    return g @ g                           # 单位权重矩阵(演示)

res = minimize(gmm_obj, x0=[0.0], method='Nelder-Mead')
print(f"估计的平均消费增长 delta0 = {res.x[0]:.4f}")
print(f"消费增长样本标准差 = {dep.std():.4f}")
# 过度识别检验: 工具个数(2个排除性工具) > 待估参数(1个)
print(f"矩条件个数 = {Z.shape[1]}, 待估参数 = 1, 自由度={Z.shape[1]-1}")

*==============================================================*
* 欧拉方程 GMM 估计 (Runkle 1991) — Stata 版
* 数据: 家庭面板 id-year, 变量 lnc=log消费, lny=log收入
*==============================================================*
clear all
set more off
use "panel_lifecycle.dta", clear
xtset id year

* 消费增长
gen dlnc = D.lnc
gen dlny = D.lny

* 欧拉方程: dlnc_{t+1} = a + e
* 工具: 滞后2期的消费/收入增长 (t期信息集, 与e_{t+1}正交)
* 用滞后变量作工具, 规避当期测量误差与预期误差的同期相关
ivregress gmm dlnc L2.dlnc L2.dlny, wmatrix(robust)
estimates store euler_iv

* 过度识别 Hansen J 检验 ( estat overid )
estat overid

* 若要识别 sigma: 结合已知 r 与校准 beta,
* 由 delta0 = (1/sigma)*(ln beta + ln(1+r)) 反推 sigma = (ln beta + ln(1+r))/delta0