保险需求与风险偏好的结构估计：从 CRRA 到逆向选择

"""
代码块1：保险需求模型 + 混合Logit估计
  U_j = (1-theta)u(w-p_j) + theta u(w-p_j-(1-q_j)L) + eps_j
  CRRA u(c)=c^(1-gamma)/(1-gamma)；theta 与 gamma 均随人异质
"""
import numpy as np
from scipy.optimize import minimize
from scipy.stats import norm
np.random.seed(20240912)

N, w, L = 2000, 10.0, 1.0
p0, q0 = 0.0, 0.0
p1, q1 = 2.0, 0.8                       # 保险合同：保费2，赔付80%

# 真实风险类型 theta ~ Logit-Normal（均值~0.25）
logit_t = norm.ppf(0.25) + 0.8*np.random.randn(N)
theta = 1/(1+np.exp(-logit_t))
# 真实风险厌恶 gamma ~ N(1.5, 0.5^2)
gamma = np.clip(1.5 + 0.5*np.random.randn(N), 0.1, 5)

def crra(c, g):
    c = np.clip(c, 1e-8, None)
    return np.where(np.abs(g-1)<1e-8, np.log(c), c**(1-g)/(1-g))

def EU(theta_i, g, p, q):
    return (1-theta_i)*crra(w-p, g) + theta_i*crra(w-p-(1-q)*L, g)

EU0 = EU(theta, gamma, p0, q0)
EU1 = EU(theta, gamma, p1, q1)
P_buy = 1/(1+np.exp(-0.5*(EU1-EU0)))     # 尺度=0.5
buy = (np.random.rand(N) < P_buy).astype(int)
print("总体购买率 =", round(buy.mean(),3))

# 混合Logit估计：在R个gamma draws上积分
def neg_ll(params, t_obs, choice):
    g_mean, g_sd = params
    if g_sd <= 1e-3: return 1e10
    lam, R = 0.5, 30
    draws = np.clip(g_mean + g_sd*np.random.randn(R,N), 0.05, 5)
    P = np.zeros(N)
    for r in range(R):
        d0 = EU(t_obs, draws[r], p0, q0)
        d1 = EU(t_obs, draws[r], p1, q1)
        P += (1/(1+np.exp(-lam*(d1-d0))))/R
    P = np.clip(P, 1e-6, 1-1e-6)
    return -np.sum(choice*np.log(P)+(1-choice)*np.log(1-P))

res = minimize(neg_ll, [1.0,0.3], args=(theta,buy),
               method='Nelder-Mead', options={'maxiter':500})
print("估计 gamma 均值 =", round(res.x[0],3), "(真实 1.5)")
print("估计 gamma 标准差 =", round(res.x[1],3), "(真实 0.5)")

"""
代码块2：从 WTP 反推风险偏好 gamma（确定性等价/风险溢价法）
  WTP_i = 为全额保险愿付最高保费；u(w-WTP) = EU_no(gamma)
  已知 theta，数值求根解 gamma
"""
import numpy as np
from scipy.optimize import brentq
np.random.seed(20240912)
w, L, N = 10.0, 2.0, 5000
gamma_true = np.random.lognormal(np.log(1.5), 0.3, N)
theta = np.full(N, 0.3)

def crra(c, g):
    c = np.clip(c, 1e-8, None)
    return c**(1-g)/(1-g)

def max_premium(t, g):
    EU_no = (1-t)*crra(w,g) + t*crra(w-L,g)
    return brentq(lambda p: crra(w-p,g)-EU_no, 1e-6, w-1e-3)

WTP = np.array([max_premium(theta[i], gamma_true[i]) for i in range(N)])

def recover_gamma(wtp, t):
    def eq(g):
        if g < 0.01: return 1e6
        return crra(w-wtp,g) - ((1-t)*crra(w,g)+t*crra(w-L,g))
    return brentq(eq, 0.01, 8.0)

g_hat = np.array([recover_gamma(WTP[i], theta[i]) for i in range(N)])
print("E[L] =", theta[0]*L)
print("真实 gamma 均值 =", round(gamma_true.mean(),3))
print("反推 gamma 均值 =", round(np.nanmean(g_hat),3))

"""
代码块3：逆向选择均衡 + 社会最优反事实（Rothschild-Stiglitz）
  两类型：theta_H=0.45, theta_L=0.15；高风险占比 phi=0.35
  市场保费 p = 购买者平均风险*L（出清）；社会最优 p=人群平均风险*L
"""
import numpy as np
from scipy.optimize import brentq
np.random.seed(20240912)
w, L, gamma = 10.0, 1.0, 2.0
theta_L, theta_H, frac_H = 0.15, 0.45, 0.35
N = 10000
type_H = np.random.rand(N) < frac_H
theta_i = np.where(type_H, theta_H, theta_L)

def crra(c, g):
    c = np.clip(c, 1e-8, None); return c**(1-g)/(1-g)

def buys(p, t):
    return crra(w-p, gamma) >= (1-t)*crra(w,gamma)+t*crra(w-L,gamma)

# 市场均衡固定点
p = 0.3
for _ in range(200):
    b = buys(p, theta_i)
    if b.sum()==0: p=np.nan; break
    p_new = theta_i[b].mean()*L
    if abs(p_new-p)<1e-6: p=p_new; break
    p = 0.5*p+0.5*p_new
print("市场均衡保费 p* =", round(p,3))
b_eq = buys(p, theta_i)
print("低风险购买率 =", round(b_eq[~type_H].mean(),3))
print("高风险购买率 =", round(b_eq[type_H].mean(),3))

# 社会最优（强制按人群平均风险）
p_social = (frac_H*theta_H+(1-frac_H)*theta_L)*L
print("社会最优保费 =", round(p_social,3))
print("解读：市场保费被高风险推高，低风险退出 -> 保障不足 -> 福利损失。")
print("      强制/补贴保险可恢复覆盖、消除该福利损失。")