导入数据¶
In [1]:
import pandas as pd
c0m0full = pd.read_csv(r'D:\论文\最后一波一鼓作气\数据\c0m0full.csv')
c1m1full = pd.read_csv(r'D:\论文\最后一波一鼓作气\数据\c1m1full.csv')
hsf15full = pd.read_csv(r'D:\论文\最后一波一鼓作气\数据\hsf15full.csv')
hsf18full = pd.read_csv(r'D:\论文\最后一波一鼓作气\数据\hsf18full.csv')
city_gender_age_premium_ratio_district15 = pd.read_csv(r'D:\论文\最后一波一鼓作气\数据\city_gender_age_premium_ratio_district15.csv')
In [2]:
#删除城市样本
c1m1full = c1m1full[c1m1full['urban_nbs'] != 'Urban']
#确保15年未整合，18年整合了
c1m1full = c1m1full[(c1m1full['policyintergration2015']==0.0) & (c1m1full['policyintergration2018']==1.0)]
c1m1= c1m1full[['ID', 'c1','m1']]
#删除城市样本
c0m0full = c0m0full[c0m0full['urban_nbs'] != 'Urban']
#确保15年未整合，18年整合了
c0m0full = c0m0full[(c0m0full['policyintergration2015']==0.0) & (c0m0full['policyintergration2018']==1.0)]
c0m0= c0m0full[['ID', 'c0','m0']]
#删除城市样本
hsf15full = hsf15full[hsf15full['urban_nbs'] != 'Urban']
#确保15年未整合，18年整合了
hsf15full = hsf15full[(hsf15full['policyintergration2015']==0.0) & (hsf15full['policyintergration2018']==1.0)]
hsf15= hsf15full[['ID', 'hsf15']]
#删除城市样本
hsf18full = hsf18full[hsf18full['urban_nbs'] != 'Urban']
#确保15年未整合，18年整合了
hsf18full = hsf18full[(hsf18full['policyintergration2015']==0.0) & (hsf18full['policyintergration2018']==1.0)]
hsf18= hsf18full[['ID', 'hsf18']]
最优化方法¶
consumption-based：合并数据¶
In [3]:
data = pd.merge(c0m0, c1m1full,on="ID", how="inner")
data
Out[3]:
| ID | c0 | m0 | householdID | communityID | c1 | m1 | gender | age | marriage | ... | premium2018 | r0 | r1 | r0adjust | r1adjust | policyintergration2015 | policyintergration2018 | district | GDPgrowthrate | urban_nbs | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 64033321002 | 112.216 | 60.0 | 640333210 | 640333 | 123.670 | 0.0 | 0.0 | 59.0 | 1.0 | ... | 180.0 | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural |
| 1 | 64033327002 | 1029.200 | 6000.0 | 640333270 | 640333 | 1054.764 | 21000.0 | 0.0 | 62.0 | 1.0 | ... | 180.0 | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural |
| 2 | 64033325001 | 1062.400 | 3000.0 | 640333250 | 640333 | 78.020 | 100.0 | 1.0 | 66.0 | 1.0 | ... | 180.0 | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural |
| 3 | 64033322001 | 592.620 | 2000.0 | 640333220 | 640333 | 859.050 | 17050.0 | 1.0 | 63.0 | 1.0 | ... | 180.0 | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural |
| 4 | 64033330002 | 2058.400 | 2000.0 | 640333300 | 640333 | 4025.500 | 2000.0 | 1.0 | 59.0 | 1.0 | ... | 180.0 | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 3882 | 89676104001 | 2315.700 | 840.0 | 896761040 | 896761 | 891.420 | 8000.0 | 1.0 | 61.0 | 1.0 | ... | 220.0 | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural |
| 3883 | 89676114002 | 1935.145 | 300.0 | 896761140 | 896761 | 49.800 | 0.0 | 0.0 | 56.0 | 0.0 | ... | 220.0 | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural |
| 3884 | 89676118001 | 2466.096 | 1000.0 | 896761180 | 896761 | 661.095 | 3000.0 | 0.0 | 73.0 | 0.0 | ... | 220.0 | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural |
| 3885 | 89676115001 | 10721.940 | 800.0 | 896761150 | 896761 | 11638.260 | 2000.0 | 1.0 | 55.0 | 1.0 | ... | 220.0 | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural |
| 3886 | 89676124001 | 268.422 | 500.0 | 896761240 | 896761 | 313.242 | 101.0 | 1.0 | 69.0 | 0.0 | ... | 220.0 | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural |
3887 rows × 26 columns
In [4]:
#最优化方法——消费计算
import pandas as pd
import numpy as np
# 计算 E(1/c0) 和 E(1/c1)
E_c0_inv2 = (1/data['c0']).mean()
E_c1_inv2 = (1/data['c1']).mean()
# 计算协方差
cov_c0_inv2 = np.cov((1/data['c0']) / E_c0_inv2, (data['r0'] - data['r1']) * data['m0'] + data['premium2015'] - data['premium2018'])[0, 1]
cov_c1_inv2 = np.cov((1/data['c1']) / E_c1_inv2, (data['r0'] - data['r1']) * data['m1'] + data['premium2015'] - data['premium2018'])[0, 1]
gamma712=abs(city_gender_age_premium_ratio_district15['premium2015'].mean() - city_gender_age_premium_ratio_district15['premium2018'].mean()) + abs(0.5 * (city_gender_age_premium_ratio_district15['r0'].mean() - city_gender_age_premium_ratio_district15['r1'].mean()) * (c0m0['m0'].mean() + c1m1['m1'].mean())) + 0.5 * cov_c0_inv2 + 0.5 * cov_c1_inv2
float(gamma712)
Out[4]:
938.0511031413579
In [5]:
#异质性男性 混合截面
#最优化方法——消费计算
datamale=data[data['gender'] == 1]
import pandas as pd
import numpy as np
E_c0_inv2 = (1/datamale['c0']).mean()
E_c1_inv2 = (1/datamale['c1']).mean()
# 计算协方差
cov_c0_inv2 = np.cov(1/datamale['c0']/ E_c0_inv2, (datamale['r0'] - datamale['r1']) * datamale['m0'] + datamale['premium2015'] - datamale['premium2018'])[0, 1]
cov_c1_inv2 = np.cov(1/datamale['c1']/ E_c1_inv2, (datamale['r0'] - datamale['r1']) * datamale['m1'] + datamale['premium2015'] - datamale['premium2018'])[0, 1]
gamma722=abs(city_gender_age_premium_ratio_district15[city_gender_age_premium_ratio_district15['gender'] == 1]['premium2015'].mean() - city_gender_age_premium_ratio_district15[city_gender_age_premium_ratio_district15['gender'] == 1]['premium2018'].mean()) + abs(0.5 * (city_gender_age_premium_ratio_district15[city_gender_age_premium_ratio_district15['gender'] == 1]['r0'].mean() - city_gender_age_premium_ratio_district15[city_gender_age_premium_ratio_district15['gender'] == 1]['r1'].mean()) * (c0m0full[c0m0full['gender'] == 1]['m0'].mean() + c1m1full[c1m1full['gender'] == 1]['m1'].mean())) + 0.5 * cov_c0_inv2 + 0.5 * cov_c1_inv2
float(gamma722)
Out[5]:
899.4616551269398
In [6]:
#异质性女性 混合截面
#最优化方法——消费计算
datafemale=data[data['gender'] == 0]
import pandas as pd
import numpy as np
E_c0_inv2 = (1/datafemale['c0']).mean()
E_c1_inv2 = (1/datafemale['c1']).mean()
# 计算协方差
cov_c0_inv2 = np.cov(1/datafemale['c0'] / E_c0_inv2, (datafemale['r0'] - datafemale['r1']) * datafemale['m0'] + datafemale['premium2015'] - datafemale['premium2018'])[0, 1]
cov_c1_inv2 = np.cov(1/datafemale['c1'] / E_c1_inv2, (datafemale['r0'] - datafemale['r1']) * datafemale['m1'] + datafemale['premium2015'] - datafemale['premium2018'])[0, 1]
gamma732=abs(city_gender_age_premium_ratio_district15[city_gender_age_premium_ratio_district15['gender'] == 0]['premium2015'].mean() - city_gender_age_premium_ratio_district15[city_gender_age_premium_ratio_district15['gender'] == 0]['premium2018'].mean()) + abs(0.5 * (city_gender_age_premium_ratio_district15[city_gender_age_premium_ratio_district15['gender'] == 0]['r0'].mean() - city_gender_age_premium_ratio_district15[city_gender_age_premium_ratio_district15['gender'] == 0]['r1'].mean()) * (c0m0full[c0m0full['gender'] == 0]['m0'].mean() + c1m1full[c1m1full['gender'] == 0]['m1'].mean())) + 0.5 * cov_c0_inv2 + 0.5 * cov_c1_inv2
float(gamma732)
Out[6]:
977.47756510789
In [7]:
import numpy as np
import pandas as pd
# ---------- 小工具 ----------
def safe_filter(df, cond_fn):
"""对 df 应用条件；若缺列/异常则返回原 df（不筛选）"""
try:
m = cond_fn(df)
if isinstance(m, pd.Series) and len(m) == len(df):
return df.loc[m]
except Exception:
pass
return df
def safe_mean(df, col):
"""数值均值（忽略缺失）；若列不存在返回 NaN"""
if col not in df.columns:
return np.nan
return pd.to_numeric(df[col], errors="coerce").mean(skipna=True)
def safe_cov(x, y):
"""样本协方差（忽略 NaN/Inf；样本<2 返回 0）"""
z = pd.concat([x, y], axis=1).replace([np.inf, -np.inf], np.nan).dropna()
if len(z) >= 2:
return float(np.cov(z.iloc[:, 0], z.iloc[:, 1], ddof=1)[0, 1])
return 0.0
def compute_gamma_inverse(data_sub, city_sub, m0_sub, m1_sub):
"""本方法：Ec0=E[1/c0]、Ec1=E[1/c1]，协方差用 (1/c)/Ec 标准化"""
# 转数值
c0 = pd.to_numeric(data_sub["c0"], errors="coerce")
c1 = pd.to_numeric(data_sub["c1"], errors="coerce")
r0 = pd.to_numeric(data_sub["r0"], errors="coerce")
r1 = pd.to_numeric(data_sub["r1"], errors="coerce")
m0 = pd.to_numeric(data_sub["m0"], errors="coerce")
m1 = pd.to_numeric(data_sub["m1"], errors="coerce")
p15 = pd.to_numeric(data_sub["premium2015"], errors="coerce")
p18 = pd.to_numeric(data_sub["premium2018"], errors="coerce")
# 避免除以 0：先把 0 设为 NaN
inv_c0 = 1.0 / c0.replace(0, np.nan)
inv_c1 = 1.0 / c1.replace(0, np.nan)
Ec0 = inv_c0.mean()
Ec1 = inv_c1.mean()
cov0 = cov1 = 0.0
if pd.notna(Ec0) and Ec0 != 0:
x0 = inv_c0 / Ec0
y0 = (r0 - r1) * m0 + (p15 - p18)
cov0 = safe_cov(x0, y0)
if pd.notna(Ec1) and Ec1 != 0:
x1 = inv_c1 / Ec1
y1 = (r0 - r1) * m1 + (p15 - p18)
cov1 = safe_cov(x1, y1)
delta_premium = safe_mean(city_sub, "premium2015") - safe_mean(city_sub, "premium2018")
delta_r = safe_mean(city_sub, "r0") - safe_mean(city_sub, "r1")
avg_m = safe_mean(m0_sub, "m0") + safe_mean(m1_sub, "m1")
gamma = abs(delta_premium) + abs(0.5 * delta_r * avg_m) + 0.5 * cov0 + 0.5 * cov1
return float(gamma)
# ---------- 异质性条件（编号 2~22） ----------
conds = {
2: lambda d: d["gender"].eq(1), # 男
3: lambda d: d["gender"].eq(0), # 女
4: lambda d: d["marriage"].eq(1), # marriage=1
5: lambda d: d["marriage"].eq(0), # marriage=0
6: lambda d: d["kids15"].eq(1), # kids15=1
7: lambda d: d["kids15"].eq(0), # kids15=0
8: lambda d: d["age"] < 59, # age<59
9: lambda d: d["age"].between(60, 79, inclusive="both"), # 60~79
10: lambda d: d["age"] >= 80, # 80+
11: lambda d: d["district"].astype(str).str.lower().eq("east"), # east
12: lambda d: d["district"].astype(str).str.lower().eq("middle"), # middle
13: lambda d: d["district"].astype(str).str.lower().eq("west"), # west
14: lambda d: d["hsf15"] > 40, # hsf15>40
15: lambda d: d["hsf15"].between(25, 40, inclusive="both"), # 25~40
16: lambda d: d["hsf15"] < 25, # <25
17: lambda d: d["ic15"] > 35000, # ic15>35000
18: lambda d: d["ic15"].between(5000, 35000, inclusive="both"), # 5000~35000
19: lambda d: d["ic15"] < 5000, # <5000
20: lambda d: d["educationrevised"].isin([6,7,8,9,10,11]), # 教育 6-11
21: lambda d: d["educationrevised"].eq(5), # 教育 5
22: lambda d: d["educationrevised"].isin([1,2,3,4]), # 教育 1-4
}
# ---------- 批量计算：gamma622 … gamma6222 ----------
_results_inv = {}
for idx, cond_fn in conds.items():
data_sub = safe_filter(data, cond_fn)
city_sub = safe_filter(city_gender_age_premium_ratio_district15, cond_fn)
m0_sub = safe_filter(c0m0full, cond_fn)
m1_sub = safe_filter(c1m1full, cond_fn)
name = f"gamma7{idx}2" # 末尾 2 = 本“1/c”方法
_results_inv[name] = compute_gamma_inverse(data_sub, city_sub, m0_sub, m1_sub)
# 可选：注册为同名变量
globals().update(_results_inv)
# 打印核对
for idx in range(2, 23):
key = f"gamma7{idx}2"
print(f"{key} = {_results_inv.get(key, np.nan)}")
gamma722 = 899.4616551269398 gamma732 = 977.47756510789 gamma742 = 972.9368025379024 gamma752 = 632.0173940199315 gamma762 = 922.7521802184212 gamma772 = 443.37617380463996 gamma782 = 819.7923011456289 gamma792 = 1002.1821354527082 gamma7102 = 865.0747449209646 gamma7112 = 957.6806214654343 gamma7122 = 988.9734901867552 gamma7132 = 795.801408446174 gamma7142 = 859.7196389259227 gamma7152 = 1055.0313617071943 gamma7162 = 832.3105919566708 gamma7172 = 848.0542634313481 gamma7182 = 708.1596131038498 gamma7192 = 1031.3260150401989 gamma7202 = 891.426034633548 gamma7212 = 1027.203575599706 gamma7222 = 905.2606829307243
health-based：合并数据¶
In [8]:
dataa = pd.merge(c0m0, c1m1full,on="ID", how="inner")
dataa
Out[8]:
| ID | c0 | m0 | householdID | communityID | c1 | m1 | gender | age | marriage | ... | premium2018 | r0 | r1 | r0adjust | r1adjust | policyintergration2015 | policyintergration2018 | district | GDPgrowthrate | urban_nbs | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 64033321002 | 112.216 | 60.0 | 640333210 | 640333 | 123.670 | 0.0 | 0.0 | 59.0 | 1.0 | ... | 180.0 | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural |
| 1 | 64033327002 | 1029.200 | 6000.0 | 640333270 | 640333 | 1054.764 | 21000.0 | 0.0 | 62.0 | 1.0 | ... | 180.0 | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural |
| 2 | 64033325001 | 1062.400 | 3000.0 | 640333250 | 640333 | 78.020 | 100.0 | 1.0 | 66.0 | 1.0 | ... | 180.0 | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural |
| 3 | 64033322001 | 592.620 | 2000.0 | 640333220 | 640333 | 859.050 | 17050.0 | 1.0 | 63.0 | 1.0 | ... | 180.0 | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural |
| 4 | 64033330002 | 2058.400 | 2000.0 | 640333300 | 640333 | 4025.500 | 2000.0 | 1.0 | 59.0 | 1.0 | ... | 180.0 | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 3882 | 89676104001 | 2315.700 | 840.0 | 896761040 | 896761 | 891.420 | 8000.0 | 1.0 | 61.0 | 1.0 | ... | 220.0 | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural |
| 3883 | 89676114002 | 1935.145 | 300.0 | 896761140 | 896761 | 49.800 | 0.0 | 0.0 | 56.0 | 0.0 | ... | 220.0 | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural |
| 3884 | 89676118001 | 2466.096 | 1000.0 | 896761180 | 896761 | 661.095 | 3000.0 | 0.0 | 73.0 | 0.0 | ... | 220.0 | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural |
| 3885 | 89676115001 | 10721.940 | 800.0 | 896761150 | 896761 | 11638.260 | 2000.0 | 1.0 | 55.0 | 1.0 | ... | 220.0 | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural |
| 3886 | 89676124001 | 268.422 | 500.0 | 896761240 | 896761 | 313.242 | 101.0 | 1.0 | 69.0 | 0.0 | ... | 220.0 | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural |
3887 rows × 26 columns
In [9]:
#计算dh/dm 15
import pandas as pd
import statsmodels.api as sm
e1 = pd.merge(c0m0,hsf15,on="ID",how="inner")
e1= e1[['m0','hsf15']].copy()
# 删除包含 NaN 或 inf 的行
e1= e1.replace([np.inf, -np.inf], np.nan).dropna()
# 删除包含 0 的行
e1 = e1[(e1['hsf15']!=0) & (e1['m0']!=0)]
# 自变量（X）和因变量（Y）
X = e1['m0']
Y = e1['hsf15']
# 在 X 中添加常数项，以便进行 OLS 回归
X = sm.add_constant(X)
# 拟合 OLS 回归模型
model = sm.OLS(Y, X).fit()
# 输出回归结果
print(model.summary())
# 提取回归系数
coefficients = model.params
# 保存特定自变量的回归系数
h_m15 = coefficients['m0']
h_m15
OLS Regression Results
==============================================================================
Dep. Variable: hsf15 R-squared: 0.000
Model: OLS Adj. R-squared: 0.000
Method: Least Squares F-statistic: 1.297
Date: Mon, 13 Oct 2025 Prob (F-statistic): 0.255
Time: 23:05:55 Log-Likelihood: -11963.
No. Observations: 3113 AIC: 2.393e+04
Df Residuals: 3111 BIC: 2.394e+04
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 33.6007 0.216 155.587 0.000 33.177 34.024
m0 -1.706e-05 1.5e-05 -1.139 0.255 -4.64e-05 1.23e-05
==============================================================================
Omnibus: 41.096 Durbin-Watson: 1.849
Prob(Omnibus): 0.000 Jarque-Bera (JB): 32.348
Skew: -0.167 Prob(JB): 9.46e-08
Kurtosis: 2.629 Cond. No. 1.54e+04
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.54e+04. This might indicate that there are
strong multicollinearity or other numerical problems.
Out[9]:
np.float64(-1.7062284276591698e-05)
In [10]:
#计算dh/dm 18
import pandas as pd
import statsmodels.api as sm
f1 = pd.merge(c1m1,hsf18,on="ID",how="inner")
f1= f1[['m1','hsf18']].copy()
# 删除包含 NaN 或 inf 的行
f1= f1.replace([np.inf, -np.inf], np.nan).dropna()
# 删除包含 0 的行
f1 = f1[(f1['hsf18']!=0) & (f1['m1']!=0)]
# 自变量（X）和因变量（Y）
X = f1['m1']
Y = f1['hsf18']
# 在 X 中添加常数项，以便进行 OLS 回归
X = sm.add_constant(X)
# 拟合 OLS 回归模型
model = sm.OLS(Y, X).fit()
# 输出回归结果
print(model.summary())
# 提取回归系数
coefficients = model.params
# 保存特定自变量的回归系数
h_m18 = coefficients['m1']
h_m18
OLS Regression Results
==============================================================================
Dep. Variable: hsf18 R-squared: 0.000
Model: OLS Adj. R-squared: 0.000
Method: Least Squares F-statistic: 1.579
Date: Mon, 13 Oct 2025 Prob (F-statistic): 0.209
Time: 23:05:55 Log-Likelihood: -26371.
No. Observations: 5630 AIC: 5.275e+04
Df Residuals: 5628 BIC: 5.276e+04
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 51.5258 0.359 143.352 0.000 50.821 52.230
m1 -1.359e-05 1.08e-05 -1.257 0.209 -3.48e-05 7.61e-06
==============================================================================
Omnibus: 14923.644 Durbin-Watson: 1.748
Prob(Omnibus): 0.000 Jarque-Bera (JB): 441.401
Skew: -0.235 Prob(JB): 1.42e-96
Kurtosis: 1.711 Cond. No. 3.42e+04
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 3.42e+04. This might indicate that there are
strong multicollinearity or other numerical problems.
Out[10]:
np.float64(-1.3590936773055398e-05)
In [11]:
#最优化方法——健康计算
import pandas as pd
import numpy as np
# 计算 E(1/c0) 和 E(1/c1)
E_c0_inv2 = (1/dataa['c0']).mean()
E_c1_inv2 = (1/dataa['c1']).mean()
# 计算协方差
cov_c0_inv2 = np.cov((0.019743 * h_m15)/ (E_c0_inv2 * dataa['r0']), (dataa['r0'] - dataa['r1']) * dataa['m0'] + dataa['premium2015'] - dataa['premium2018'])[0, 1]
cov_c1_inv2 = np.cov((0.019743 * h_m18)/ (E_c1_inv2 * dataa['r1']), (dataa['r0'] - dataa['r1']) * dataa['m1'] + dataa['premium2015'] - dataa['premium2018'])[0, 1]
gamma713=abs(city_gender_age_premium_ratio_district15['premium2015'].mean() - city_gender_age_premium_ratio_district15['premium2018'].mean()) + abs(0.5 * (city_gender_age_premium_ratio_district15['r0'].mean() - city_gender_age_premium_ratio_district15['r1'].mean()) * (c0m0['m0'].mean() + c1m1['m1'].mean())) + 0.5 * cov_c0_inv2 + 0.5 * cov_c1_inv2
float(gamma713)
Out[11]:
761.8113887118584
In [12]:
#表6健康based异质性—男性
dataamale=dataa[dataa['gender'] == 1]
import pandas as pd
import numpy as np
E_c0_inv2 = (1/dataamale['c0']).mean()
E_c1_inv2 = (1/dataamale['c1']).mean()
# 计算协方差
cov_c0_inv2 = np.cov((0.019743 * h_m15)/ (E_c0_inv2 * dataamale['r0']), (dataamale['r0'] - dataamale['r1']) * dataamale['m0'] + dataamale['premium2015'] - dataamale['premium2018'])[0, 1]
cov_c1_inv2 = np.cov((0.019743 * h_m18)/ (E_c1_inv2 * dataamale['r1']), (dataamale['r0'] - dataamale['r1']) * dataamale['m1'] + dataamale['premium2015'] - dataamale['premium2018'])[0, 1]
#最优化方法——健康计算（表6health based混合截面）
gamma723=abs(city_gender_age_premium_ratio_district15[city_gender_age_premium_ratio_district15['gender'] == 1]['premium2015'].mean() - city_gender_age_premium_ratio_district15[city_gender_age_premium_ratio_district15['gender'] == 1]['premium2018'].mean()) + abs(0.5 * (city_gender_age_premium_ratio_district15[city_gender_age_premium_ratio_district15['gender'] == 1]['r0'].mean() - city_gender_age_premium_ratio_district15[city_gender_age_premium_ratio_district15['gender'] == 1]['r1'].mean()) * (c0m0full[c0m0full['gender'] == 1]['m0'].mean() + c1m1full[c1m1full['gender'] == 1]['m1'].mean())) + 0.5 * cov_c0_inv2 + 0.5 * cov_c1_inv2
float(gamma723)
Out[12]:
772.2338794777837
In [13]:
#健康based异质性—女性
dataafemale=dataa[dataa['gender'] == 0]
import pandas as pd
import numpy as np
E_c0_inv2 = (1/dataafemale['c0']).mean()
E_c1_inv2 = (1/dataafemale['c1']).mean()
# 计算协方差
cov_c0_inv2 = np.cov((0.019743 * h_m15)/ (E_c0_inv2 * dataafemale['r0']), (dataafemale['r0'] - dataafemale['r1']) * dataafemale['m0'] + dataafemale['premium2015'] - dataafemale['premium2018'])[0, 1]
cov_c1_inv2 = np.cov((0.019743 * h_m18)/ (E_c1_inv2 * dataafemale['r1']), (dataafemale['r0'] - dataafemale['r1']) * dataafemale['m1'] + dataafemale['premium2015'] - dataafemale['premium2018'])[0, 1]
#最优化方法——健康计算（表6health based混合截面）
gamma723=abs(city_gender_age_premium_ratio_district15[city_gender_age_premium_ratio_district15['gender'] == 0]['premium2015'].mean() - city_gender_age_premium_ratio_district15[city_gender_age_premium_ratio_district15['gender'] == 0]['premium2018'].mean()) + abs(0.5 * (city_gender_age_premium_ratio_district15[city_gender_age_premium_ratio_district15['gender'] == 0]['r0'].mean() - city_gender_age_premium_ratio_district15[city_gender_age_premium_ratio_district15['gender'] == 0]['r1'].mean()) * (c0m0full[c0m0full['gender'] == 0]['m0'].mean() + c1m1full[c1m1full['gender'] == 0]['m1'].mean())) + 0.5 * cov_c0_inv2 + 0.5 * cov_c1_inv2
float(gamma723)
Out[13]:
750.3489588282599
In [14]:
import numpy as np
import pandas as pd
PHI = 0.019743 # phi_tilde
# ---------- 小工具 ----------
def safe_filter(df, cond_fn):
"""对 df 应用条件；若缺列/异常则返回原 df（不筛选）"""
try:
m = cond_fn(df)
if isinstance(m, pd.Series) and len(m) == len(df):
return df.loc[m]
except Exception:
pass
return df
def safe_mean(df, col):
"""数值均值（忽略缺失）；若列不存在返回 NaN"""
if col not in df.columns:
return np.nan
return pd.to_numeric(df[col], errors="coerce").mean(skipna=True)
def safe_cov(x, y):
"""样本协方差（忽略 NaN/Inf；样本<2 返回 0）"""
z = pd.concat([x, y], axis=1).replace([np.inf, -np.inf], np.nan).dropna()
if len(z) >= 2:
return float(np.cov(z.iloc[:, 0], z.iloc[:, 1], ddof=1)[0, 1])
return 0.0
def _align_health_for_subset(dataa_full, data_sub, h_full, id_col="ID"):
"""
将全样本的 h_m15/h_m18 与 data_sub 对齐：
- 若 h_full 是 Series 且索引=dataa_full.index，则按 index 对齐
- 若 h_full 索引是 ID，则按 ID 对齐
- 若 h_full 是 array/list，则包成 Series(index=dataa_full.index) 后再对齐
"""
if isinstance(h_full, pd.Series):
if h_full.index.equals(dataa_full.index):
return h_full.loc[data_sub.index]
if id_col in data_sub.columns and h_full.index.isin(data_sub[id_col]).any():
s = h_full.reindex(data_sub[id_col])
s.index = data_sub.index
return s
try:
return h_full.loc[data_sub.index]
except Exception:
return pd.Series(np.nan, index=data_sub.index)
# array-like
try:
base = pd.Series(h_full, index=dataa_full.index)
return base.loc[data_sub.index]
except Exception:
return pd.Series(np.nan, index=data_sub.index)
def compute_gamma_health_inverse(dataa_full, data_sub, city_sub, m0_sub, m1_sub, h_m15, h_m18):
"""
本方法（表6）：E[1/c] 标准化 + 健康项
Ec0 = E[1/c0], Ec1 = E[1/c1]
cov0 = Cov( (PHI * h_m15)/(Ec0 * r0), (r0-r1)*m0 + p15 - p18 )
cov1 = Cov( (PHI * h_m18)/(Ec1 * r1), (r0-r1)*m1 + p15 - p18 )
γ = |Δpremium| + |0.5*Δr*(E[m0]+E[m1])| + 0.5*(cov0 + cov1)
"""
# 数值列
c0 = pd.to_numeric(data_sub["c0"], errors="coerce")
c1 = pd.to_numeric(data_sub["c1"], errors="coerce")
r0 = pd.to_numeric(data_sub["r0"], errors="coerce")
r1 = pd.to_numeric(data_sub["r1"], errors="coerce")
m0 = pd.to_numeric(data_sub["m0"], errors="coerce")
m1 = pd.to_numeric(data_sub["m1"], errors="coerce")
p15 = pd.to_numeric(data_sub["premium2015"], errors="coerce")
p18 = pd.to_numeric(data_sub["premium2018"], errors="coerce")
# 避免除以 0：把 0 视为缺失
inv_c0 = 1.0 / c0.replace(0, np.nan)
inv_c1 = 1.0 / c1.replace(0, np.nan)
Ec0 = inv_c0.mean()
Ec1 = inv_c1.mean()
# 对齐 h_m15 / h_m18
h15 = _align_health_for_subset(dataa_full, data_sub, h_m15)
h18 = _align_health_for_subset(dataa_full, data_sub, h_m18)
cov0 = cov1 = 0.0
if pd.notna(Ec0) and Ec0 != 0:
a0 = (PHI * h15) / (Ec0 * r0) # 可能出现 inf/NaN，safe_cov 会处理
y0 = (r0 - r1) * m0 + (p15 - p18)
cov0 = safe_cov(a0, y0)
if pd.notna(Ec1) and Ec1 != 0:
a1 = (PHI * h18) / (Ec1 * r1)
y1 = (r0 - r1) * m1 + (p15 - p18)
cov1 = safe_cov(a1, y1)
delta_premium = safe_mean(city_sub, "premium2015") - safe_mean(city_sub, "premium2018")
delta_r = safe_mean(city_sub, "r0") - safe_mean(city_sub, "r1")
avg_m = safe_mean(m0_sub, "m0") + safe_mean(m1_sub, "m1")
gamma = abs(delta_premium) + abs(0.5 * delta_r * avg_m) + 0.5 * cov0 + 0.5 * cov1
return float(gamma)
# ---------- 异质性条件（2~22） ----------
conds = {
2: lambda d: d["gender"].eq(1), # 男
3: lambda d: d["gender"].eq(0), # 女
4: lambda d: d["marriage"].eq(1), # marriage=1
5: lambda d: d["marriage"].eq(0), # marriage=0
6: lambda d: d["kids15"].eq(1), # kids15=1
7: lambda d: d["kids15"].eq(0), # kids15=0
8: lambda d: d["age"] < 59, # age<59
9: lambda d: d["age"].between(60, 79, inclusive="both"), # 60~79
10: lambda d: d["age"] >= 80, # 80+
11: lambda d: d["district"].astype(str).str.lower().eq("east"), # east
12: lambda d: d["district"].astype(str).str.lower().eq("middle"), # middle
13: lambda d: d["district"].astype(str).str.lower().eq("west"), # west
14: lambda d: d["hsf15"] > 40, # hsf15>40
15: lambda d: d["hsf15"].between(25, 40, inclusive="both"), # 25~40
16: lambda d: d["hsf15"] < 25, # <25
17: lambda d: d["ic15"] > 35000, # ic15>35000
18: lambda d: d["ic15"].between(5000, 35000, inclusive="both"), # 5000~35000
19: lambda d: d["ic15"] < 5000, # <5000
20: lambda d: d["educationrevised"].isin([6,7,8,9,10,11]), # 教育 6-11
21: lambda d: d["educationrevised"].eq(5), # 教育 5
22: lambda d: d["educationrevised"].isin([1,2,3,4]), # 教育 1-4
}
# ---------- 批量计算：gamma623 … gamma6223 ----------
_results_h = {}
for idx, cond_fn in conds.items():
data_sub = safe_filter(dataa, cond_fn)
city_sub = safe_filter(city_gender_age_premium_ratio_district15, cond_fn)
m0_sub = safe_filter(c0m0full, cond_fn)
m1_sub = safe_filter(c1m1full, cond_fn)
name = f"gamma7{idx}3"
_results_h[name] = compute_gamma_health_inverse(
dataa_full=dataa, data_sub=data_sub,
city_sub=city_sub, m0_sub=m0_sub, m1_sub=m1_sub,
h_m15=h_m15, h_m18=h_m18
)
# 可选：注册为同名变量
globals().update(_results_h)
# 打印核对
for idx in range(2, 23):
key = f"gamma7{idx}3"
print(f"{key} = {_results_h.get(key, np.nan)}")
gamma723 = 772.2338794777837 gamma733 = 750.3489588282599 gamma743 = 796.9787081494184 gamma753 = 533.8862332618123 gamma763 = 745.7448176611197 gamma773 = 413.94594318410475 gamma783 = 726.4541126233465 gamma793 = 803.6234089015956 gamma7103 = 633.1529566782665 gamma7113 = 719.9252752511742 gamma7123 = 776.4728971685053 gamma7133 = 721.1091543132715 gamma7143 = 688.8018986525792 gamma7153 = 851.6498440130374 gamma7163 = 697.3818194064443 gamma7173 = 687.8451760027752 gamma7183 = 634.6722141478845 gamma7193 = 812.2048673741209 gamma7203 = 731.9214995342486 gamma7213 = 852.7165500137347 gamma7223 = 734.6915439920137
用完全信息法求解¶
In [17]:
import numpy as np
import pandas as pd
# 常量
PHI = 0.019743
c0_bar = pd.to_numeric(c0m0["c0"], errors="coerce").mean(skipna=True)
c1_bar = pd.to_numeric(c1m1["c1"], errors="coerce").mean(skipna=True)
h0_bar = pd.to_numeric(hsf15["hsf15"], errors="coerce").mean(skipna=True)
h1_bar = pd.to_numeric(hsf18["hsf18"], errors="coerce").mean(skipna=True)
gamma711 = c1_bar * (PHI * (h1_bar - h0_bar) - np.log(c0_bar / c1_bar))
print(gamma711)
699.9028412979085
In [54]:
#异质性计算男性
import numpy as np
import pandas as pd
# 参数
PHI = 0.019743
# 取各列的均值（忽略缺失）
c0_bar = pd.to_numeric(c0m0full[c0m0full['gender'] == 1]["c0"], errors="coerce").mean(skipna=True)
c1_bar = pd.to_numeric(c1m1full[c1m1full['gender'] == 1]["c1"], errors="coerce").mean(skipna=True)
h0_bar = pd.to_numeric(hsf15full[hsf15full['gender'] == 1]["hsf15"], errors="coerce").mean(skipna=True)
h1_bar = pd.to_numeric(hsf18full[hsf18full['gender'] == 1]["hsf18"], errors="coerce").mean(skipna=True)
gamma721 = c1_bar * (PHI * (h1_bar - h0_bar) - np.log(c0_bar / c1_bar))
print(gamma721)
797.625981440441
In [55]:
#异质性计算女性
import numpy as np
import pandas as pd
# 参数
PHI = 0.019743
# 取各列的均值（忽略缺失）
c0_bar = pd.to_numeric(c0m0full[c0m0full['gender'] == 0]["c0"], errors="coerce").mean(skipna=True)
c1_bar = pd.to_numeric(c1m1full[c1m1full['gender'] == 0]["c1"], errors="coerce").mean(skipna=True)
h0_bar = pd.to_numeric(hsf15full[hsf15full['gender'] == 0]["hsf15"], errors="coerce").mean(skipna=True)
h1_bar = pd.to_numeric(hsf18full[hsf18full['gender'] == 0]["hsf18"], errors="coerce").mean(skipna=True)
gamma731 = c1_bar * (PHI * (h1_bar - h0_bar) - np.log(c0_bar / c1_bar))
print(gamma731)
627.7008677242283
In [15]:
import numpy as np
import pandas as pd
PHI = 0.019743 # 可按需修改
# --- 异质性条件 ---
conds = {
2: lambda d: d["gender"].eq(1), # 男
3: lambda d: d["gender"].eq(0), # 女
4: lambda d: d["marriage"].eq(1), # marriage=1
5: lambda d: d["marriage"].eq(0), # marriage=0
6: lambda d: d["kids15"].eq(1), # kids15=1
7: lambda d: d["kids15"].eq(0), # kids15=0
8: lambda d: d["age"] < 59, # age<59
9: lambda d: d["age"].between(60, 79, inclusive="both"), # 60~79
10: lambda d: d["age"] >= 80, # 80+
11: lambda d: d["district"].astype(str).str.lower().eq("east"), # east
12: lambda d: d["district"].astype(str).str.lower().eq("middle"), # middle
13: lambda d: d["district"].astype(str).str.lower().eq("west"), # west
14: lambda d: d["hsf15"] > 40, # hsf15>40
15: lambda d: d["hsf15"].between(25, 40, inclusive="both"), # 25~40
16: lambda d: d["hsf15"] < 25, # <25
17: lambda d: d["ic15"] > 35000, # ic15>35000
18: lambda d: d["ic15"].between(5000, 35000, inclusive="both"), # 5000~35000
19: lambda d: d["ic15"] < 5000, # <5000
20: lambda d: d["educationrevised"].isin([6,7,8,9,10,11]), # 教育 6-11
21: lambda d: d["educationrevised"].eq(5), # 教育 5
22: lambda d: d["educationrevised"].isin([1,2,3,4]), # 教育 1-4
}
# --- 工具函数（缺列时跳过筛选，保证稳健） ---
def safe_filter(df: pd.DataFrame, cond_fn):
try:
mask = cond_fn(df)
if isinstance(mask, pd.Series) and len(mask) == len(df):
return df.loc[mask]
except Exception:
pass
return df
def mean_after_filter(df: pd.DataFrame, col: str, cond_fn):
if col not in df.columns:
return np.nan
sub = safe_filter(df, cond_fn)
return pd.to_numeric(sub[col], errors="coerce").mean(skipna=True)
def gamma_from_means(cond_fn):
# 四个“full”表分别取均值（按各自能支持的条件筛选）
c0_bar = mean_after_filter(c0m0full, "c0", cond_fn)
c1_bar = mean_after_filter(c1m1full, "c1", cond_fn)
h0_bar = mean_after_filter(hsf15full, "hsf15", cond_fn)
h1_bar = mean_after_filter(hsf18full, "hsf18", cond_fn)
if any(pd.isna([c0_bar, c1_bar, h0_bar, h1_bar])):
return float("nan")
if c0_bar <= 0 or c1_bar <= 0:
return float("nan")
return float(c1_bar * (PHI * (h1_bar - h0_bar) - np.log(c0_bar / c1_bar)))
# --- 批量计算：gamma621 … gamma6221 ---
_results = {}
for idx, cond_fn in conds.items():
name = f"gamma7{idx}1"
_results[name] = gamma_from_means(cond_fn)
# 可选：注册为同名全局变量
globals().update(_results)
# 打印核对
for idx in range(2, 23):
key = f"gamma7{idx}1"
print(f"{key} = {_results.get(key, np.nan)}")
gamma721 = 797.625981440441 gamma731 = 627.7008677242283 gamma741 = 692.7487107692193 gamma751 = 348.44227070849587 gamma761 = 688.8555520139241 gamma771 = -23.245473424548475 gamma781 = 1126.7781616036805 gamma791 = 358.57108011726496 gamma7101 = -187.75052410470894 gamma7111 = 1287.819851099903 gamma7121 = 747.449339068624 gamma7131 = 256.66879269599735 gamma7141 = 990.2027469224145 gamma7151 = 713.1507005708951 gamma7161 = 391.0955581304712 gamma7171 = 728.2165712077158 gamma7181 = 571.370921047391 gamma7191 = 759.5527516130173 gamma7201 = 1951.3380773108386 gamma7211 = 1045.9301253498531 gamma7221 = 489.0048551091307
In [18]:
# -*- coding: utf-8 -*-
import pandas as pd
# 1) 行索引与数据
rows = [
"全样本",
"男性","女性",
"有配偶","无配偶",
"有子女","无子女",
"小于59 岁","60 岁—79 岁","80 岁及以上",
"东部","中部","西部",
"健康状况较好","健康状况中等","健康状况较差",
"较高收入","中等收入","较低收入",
"教育程度较高","教育程度中等","教育程度较低",
]
data = [
[gamma711, gamma712, gamma713],
[gamma721, gamma722, gamma723],
[gamma731, gamma732, gamma733],
[gamma741, gamma742, gamma743],
[gamma751, gamma752, gamma753],
[gamma761, gamma762, gamma763],
[gamma771, gamma772, gamma773],
[gamma781, gamma782, gamma783],
[gamma791, gamma792, gamma793],
[gamma7101, gamma7102, gamma7103],
[gamma7111, gamma7112, gamma7113],
[gamma7121, gamma7122, gamma7123],
[gamma7131, gamma7132, gamma7133],
[gamma7141, gamma7142, gamma7143],
[gamma7151, gamma7152, gamma7153],
[gamma7161, gamma7162, gamma7163],
[gamma7171, gamma7172, gamma7173],
[gamma7181, gamma7182, gamma7183],
[gamma7191, gamma7192, gamma7193],
[gamma7201, gamma7202, gamma7203],
[gamma7211, gamma7212, gamma7213],
[gamma7221, gamma7222, gamma7223],
]
# 2) 多级列索引
cols = pd.MultiIndex.from_tuples([
("完全信息方法",""),
("最优化方法","仅假设效用函数\n的消费部分"),
("最优化方法","仅假设效用函数\n的健康部分"),
])
df = pd.DataFrame(data, index=rows, columns=cols)
# 3) 分组起始行（加粗横线）
group_starts = {
"男性", # 性别组
"有配偶", # 婚姻组
"有子女", # 子女组
"45 岁—59 岁", # 年龄组
"东部", # 地区组
"健康状况较好", # 健康组
"较高收入", # 收入组
"教育程度较高" # 教育组
}
def row_borders(row):
label = row.name
if label in group_starts:
return ['border-top: 2px solid #4a4a4a'] * len(row)
return [''] * len(row)
# 4) 样式与展示
styler = (
df.style
.set_table_styles([
{'selector': 'th.col_heading.level0',
'props': [('font-weight', '700'),
('border-bottom','1px solid #4a4a4a')]},
{'selector': 'th.col_heading.level1',
'props': [('font-weight', '700')]},
{'selector': 'th.row_heading',
'props': [('font-weight', '700')]},
{'selector': 'table',
'props': [('border-collapse','collapse'),
('font-family','-apple-system,BlinkMacSystemFont,Segoe UI,Roboto,PingFang SC,Helvetica,Arial')]}
])
.format(precision=0)
.set_properties(**{
'text-align': 'center',
'padding': '6px',
'border':'1px solid #a0a0a0'
})
.apply(row_borders, axis=1)
)
# 在 Jupyter 中显示
styler
Out[18]:
| 完全信息方法 | 最优化方法 | ||
|---|---|---|---|
| 仅假设效用函数 的消费部分 | 仅假设效用函数 的健康部分 | ||
| 全样本 | 700 | 938 | 762 |
| 男性 | 798 | 899 | 772 |
| 女性 | 628 | 977 | 750 |
| 有配偶 | 693 | 973 | 797 |
| 无配偶 | 348 | 632 | 534 |
| 有子女 | 689 | 923 | 746 |
| 无子女 | -23 | 443 | 414 |
| 小于59 岁 | 1127 | 820 | 726 |
| 60 岁—79 岁 | 359 | 1002 | 804 |
| 80 岁及以上 | -188 | 865 | 633 |
| 东部 | 1288 | 958 | 720 |
| 中部 | 747 | 989 | 776 |
| 西部 | 257 | 796 | 721 |
| 健康状况较好 | 990 | 860 | 689 |
| 健康状况中等 | 713 | 1055 | 852 |
| 健康状况较差 | 391 | 832 | 697 |
| 较高收入 | 728 | 848 | 688 |
| 中等收入 | 571 | 708 | 635 |
| 较低收入 | 760 | 1031 | 812 |
| 教育程度较高 | 1951 | 891 | 732 |
| 教育程度中等 | 1046 | 1027 | 853 |
| 教育程度较低 | 489 | 905 | 735 |