导入数据¶
In [3]:
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 [4]:
#删除城市样本
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 [5]:
data = pd.merge(c0m0, c1m1full,on="ID", how="inner")
data
Out[5]:
| 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 [6]:
#最优化方法——消费计算
import pandas as pd
import numpy as np
# 计算 E(c0^(-3)) 和 E(c1^(-3))
E_c0_inv2 = (data['c0']**(-3)).mean()
E_c1_inv2 = (data['c1']**(-3)).mean()
# 计算协方差
cov_c0_inv2 = np.cov(data['c0']**(-3) / E_c0_inv2, (data['r0'] - data['r1']) * data['m0'] + data['premium2015'] - data['premium2018'])[0, 1]
cov_c1_inv2 = np.cov(data['c1']**(-3) / E_c1_inv2, (data['r0'] - data['r1']) * data['m1'] + data['premium2015'] - data['premium2018'])[0, 1]
#最优化方法——消费计算（表2consumption based混合截面）
gamma512=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(gamma512)
Out[6]:
990.1055666990669
In [7]:
#异质性男性 混合截面
#最优化方法——消费计算
datamale=data[data['gender'] == 1]
import pandas as pd
import numpy as np
# 计算 E(c0^(-3)) 和 E(c1^(-3))
E_c0_inv2 = (datamale['c0']**(-3)).mean()
E_c1_inv2 = (datamale['c1']**(-3)).mean()
# 计算协方差
cov_c0_inv2 = np.cov(datamale['c0']**(-3) / E_c0_inv2, (datamale['r0'] - datamale['r1']) * datamale['m0'] + datamale['premium2015'] - datamale['premium2018'])[0, 1]
cov_c1_inv2 = np.cov(datamale['c1']**(-3) / E_c1_inv2, (datamale['r0'] - datamale['r1']) * datamale['m1'] + datamale['premium2015'] - datamale['premium2018'])[0, 1]
gamma522=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(gamma522)
Out[7]:
892.3622174949148
In [8]:
#异质性女性 混合截面
#最优化方法——消费计算
datafemale=data[data['gender'] == 0]
import pandas as pd
import numpy as np
# 计算 E(c0^(-3)) 和 E(c1^(-3))
E_c0_inv2 = (datafemale['c0']**(-3)).mean()
E_c1_inv2 = (datafemale['c1']**(-3)).mean()
# 计算协方差
cov_c0_inv2 = np.cov(datafemale['c0']**(-3) / E_c0_inv2, (datafemale['r0'] - datafemale['r1']) * datafemale['m0'] + datafemale['premium2015'] - datafemale['premium2018'])[0, 1]
cov_c1_inv2 = np.cov(datafemale['c1']**(-3) / E_c1_inv2, (datafemale['r0'] - datafemale['r1']) * datafemale['m1'] + datafemale['premium2015'] - datafemale['premium2018'])[0, 1]
gamma532=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(gamma532)
Out[8]:
1068.9200342165386
In [9]:
import numpy as np
import pandas as pd
# —— 安全筛选：如果某个子表缺少相应列，就不筛选它（防止 KeyError）——
def safe_filter(df, cond_fn):
try:
mask = cond_fn(df)
# 只有当返回的是与 df 等长的布尔 Series 才进行筛选
if isinstance(mask, pd.Series) and len(mask) == len(df):
return df.loc[mask]
except Exception:
pass
return df # 缺列或报错就不筛选
# —— 计算协方差（忽略 NaN；样本量<2 时返回 0）——
def safe_cov(x, y):
tmp = pd.concat([x, y], axis=1).dropna()
if len(tmp) >= 2:
return float(np.cov(tmp.iloc[:, 0], tmp.iloc[:, 1], ddof=1)[0, 1])
return 0.0
# —— 按子样本计算 gamma ——
def compute_gamma(data_sub, city_sub, m0_sub, m1_sub):
# 期望值
E_c0_inv2 = (data_sub['c0']**(-3)).mean()
E_c1_inv2 = (data_sub['c1']**(-3)).mean()
# 协方差项
cov_c0_inv2 = 0.0
cov_c1_inv2 = 0.0
if pd.notna(E_c0_inv2) and E_c0_inv2 != 0:
x0 = (data_sub['c0']**(-3)) / E_c0_inv2
y0 = (data_sub['r0'] - data_sub['r1']) * data_sub['m0'] + (data_sub['premium2015'] - data_sub['premium2018'])
cov_c0_inv2 = safe_cov(x0, y0)
if pd.notna(E_c1_inv2) and E_c1_inv2 != 0:
x1 = (data_sub['c1']**(-3)) / E_c1_inv2
y1 = (data_sub['r0'] - data_sub['r1']) * data_sub['m1'] + (data_sub['premium2015'] - data_sub['premium2018'])
cov_c1_inv2 = safe_cov(x1, y1)
# 其它两项（来自 city_* 与 m0/m1 的“全量”表，但做相同条件筛选）
delta_premium = city_sub['premium2015'].mean() - city_sub['premium2018'].mean()
delta_r = city_sub['r0'].mean() - city_sub['r1'].mean()
avg_m = m0_sub['m0'].mean() + m1_sub['m1'].mean()
gamma = abs(delta_premium) + abs(0.5 * delta_r * avg_m) + 0.5 * cov_c0_inv2 + 0.5 * cov_c1_inv2
return float(gamma)
# —— 为每个异质性条件定义筛选函数（对 data / city_* / c0m0full / c1m1full 统一使用）——
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, # <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
}
# —— 批量计算，每个条件生成一个 gamma4{idx}2 变量 ——
_results = {}
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'gamma5{idx}2'
_results[name] = compute_gamma(data_sub, city_sub, m0_sub, m1_sub)
# 可选：把结果提升为同名变量（如 gamma442、gamma452...）
globals().update(_results)
# 查看所有结果
for idx in range(2, 23):
key = f'gamma5{idx}2'
print(f'{key} = {_results.get(key, np.nan)}')
gamma522 = 892.3622174949148 gamma532 = 1068.9200342165386 gamma542 = 1054.1372398639755 gamma552 = 591.9719665517397 gamma562 = 974.2732874972605 gamma572 = 505.38817958936147 gamma582 = 718.0145948778675 gamma592 = 1060.8951186687636 gamma5102 = 1019.9283982874636 gamma5112 = 1166.5200251872288 gamma5122 = 1029.4207535732694 gamma5132 = 789.9357996642834 gamma5142 = 906.5164898328816 gamma5152 = 1062.9499268159977 gamma5162 = 913.8399731963046 gamma5172 = 945.2397341649574 gamma5182 = 678.3323074360065 gamma5192 = 1100.7690718458314 gamma5202 = 1200.1751124667826 gamma5212 = 850.1571184651061 gamma5222 = 954.8217038429507
health-based：合并数据¶
In [10]:
dataa = pd.merge(c0m0, c1m1full,on="ID", how="inner")
dataa
Out[10]:
| 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 [11]:
#计算dh/dm 15
import pandas as pd
import numpy as np
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, 29 Dec 2025 Prob (F-statistic): 0.255
Time: 17:26:15 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[11]:
np.float64(-1.7062284276591698e-05)
In [12]:
#计算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, 29 Dec 2025 Prob (F-statistic): 0.209
Time: 17:26:18 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[12]:
np.float64(-1.3590936773055398e-05)
In [13]:
import pandas as pd
import numpy as np
# 计算 E(c0^(-3)) 和 E(c1^(-3))
E_c0_inv2 = (dataa['c0']**(-3)).mean()
E_c1_inv2 = (dataa['c1']**(-3)).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]
gamma513=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(gamma513)
Out[13]:
785.5178011684857
In [16]:
#健康based异质性探索——男性
dataamale=dataa[dataa['gender'] == 1]
import pandas as pd
import numpy as np
# 计算 E(c0^(-3)) 和 E(c1^(-3))
E_c0_inv2 = (dataamale['c0']**(-3)).mean()
E_c1_inv2 = (dataamale['c1']**(-3)).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]
#最优化方法——健康计算（表2health based混合截面）
gamma523=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(gamma523)
Out[16]:
798.2806824726699
In [15]:
import numpy as np
import pandas as pd
# ================== 1) 工具函数 ==================
phi_tilde = 0.019743 # 俺的 Φ~
def safe_filter(df, cond_fn):
"""对 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_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(data_sub, h_full, *, id_col="ID"):
"""
将全样本的 h_full（h_m15/h_m18）对齐到 data_sub 的行顺序。
支持：
- h_full 为 Series（索引=ID 或 =dataa.index）
- h_full 为 numpy 数组（长度与 dataa 相同）
"""
# 情况1：Series
if isinstance(h_full, pd.Series):
# 索引若与 dataa.index 对齐
if h_full.index.equals(dataa.index):
return h_full.loc[data_sub.index]
# 索引若是 ID（字符串/数字），用 ID 对齐
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)
# 情况2：numpy 数组或其他可迭代
try:
base = pd.Series(h_full, index=dataa.index)
return base.loc[data_sub.index]
except Exception:
return pd.Series(np.nan, index=data_sub.index)
def compute_gamma_health(data_sub, city_sub, m0_sub, m1_sub, h15_full, h18_full, phi=phi_tilde):
"""
Health-based γ：
cov0 = Cov( (φ*h_m15)/(E[c0^{-2}]*r0), (r0-r1)*m0 + premium2015 - premium2018 )
cov1 = Cov( (φ*h_m18)/(E[c1^{-2}]*r1), (r0-r1)*m1 + premium2015 - premium2018 )
γ = |Δpremium| + |0.5*Δr*(E[m0]+E[m1])| + 0.5*cov0 + 0.5*cov1
"""
# 期望（用于标准化）
E_c0_inv2 = (pd.to_numeric(data_sub["c0"], errors="coerce")**(-3)).mean()
E_c1_inv2 = (pd.to_numeric(data_sub["c1"], errors="coerce")**(-3)).mean()
# 取子样本对应的 h_m15/h_m18
h15 = _align_health_for_subset(data_sub, h15_full)
h18 = _align_health_for_subset(data_sub, h18_full)
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")
cov0 = cov1 = 0.0
if pd.notna(E_c0_inv2) and E_c0_inv2 != 0:
a0 = (phi * h15) / (E_c0_inv2 * r0)
y0 = (r0 - r1) * m0 + (p15 - p18)
cov0 = safe_cov(a0, y0)
if pd.notna(E_c1_inv2) and E_c1_inv2 != 0:
a1 = (phi * h18) / (E_c1_inv2 * r1)
y1 = (r0 - r1) * m1 + (p15 - p18)
cov1 = safe_cov(a1, y1)
delta_premium = city_sub["premium2015"].mean() - city_sub["premium2018"].mean()
delta_r = city_sub["r0"].mean() - city_sub["r1"].mean()
avg_m = m0_sub["m0"].mean() + m1_sub["m1"].mean()
gamma = abs(delta_premium) + abs(0.5 * delta_r * avg_m) + 0.5*cov0 + 0.5*cov1
return float(gamma)
# ================== 2) 条件 ==================
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, # <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
}
# ================== 3) 批量计算：gamma523 … gamma5223 ==================
_results_h = {}
for idx, cond_fn in conds.items():
d_sub = safe_filter(dataa, cond_fn)
city_s = safe_filter(city_gender_age_premium_ratio_district15, cond_fn)
m0_s = safe_filter(c0m0full, cond_fn)
m1_s = safe_filter(c1m1full, cond_fn)
name = f"gamma5{idx}3" # 末尾 3 = health-based
_results_h[name] = compute_gamma_health(d_sub, city_s, m0_s, m1_s, h_m15, h_m18, phi=phi_tilde)
# 可选：提升为同名变量
globals().update(_results_h)
# 查看结果
for idx in range(2, 23):
key = f"gamma5{idx}3"
print(f"{key} = {_results_h.get(key, np.nan)}")
gamma523 = 798.2806824726699 gamma533 = 772.6549904538863 gamma543 = 827.3541030200007 gamma553 = 548.9890006240952 gamma563 = 770.138292621716 gamma573 = 405.6688613773065 gamma583 = 750.5854229039102 gamma593 = 817.8256634765235 gamma5103 = 665.220338543296 gamma5113 = 725.9821356930689 gamma5123 = 798.5608830865662 gamma5133 = 738.729389564704 gamma5143 = 693.9697202500064 gamma5153 = 883.6557425223073 gamma5163 = 723.8315443148991 gamma5173 = 792.9577084907053 gamma5183 = 649.4122164122365 gamma5193 = 835.9700091717721 gamma5203 = 842.1899642461668 gamma5213 = 917.4985634120095 gamma5223 = 753.7055032434272
用完全信息法求解¶
In [17]:
import numpy as np
import pandas as pd
# 参数
sigmamale = 3.0
phi_tilde = 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)
B_bar = (c0_bar**(1 - sigmamale)) + (1 - sigmamale) * phi_tilde * (h0_bar - h1_bar)
cons1_bar = B_bar**(1 / (1 - sigmamale))
gamma511 = c1_bar - cons1_bar
print(gamma511)
1957.6372848184362
In [18]:
#异质性计算男性
import numpy as np
import pandas as pd
# 参数
sigmafemale = 3.0
phi_tilde = 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)
B_bar = (c0_bar**(1 - sigmafemale)) + (1 - sigmafemale) * phi_tilde * (h0_bar - h1_bar)
cons1_bar = B_bar**(1 / (1 - sigmafemale))
gamma521 = c1_bar - cons1_bar
print(gamma521)
1947.8022115045392
In [19]:
import numpy as np
import pandas as pd
# ========== 参数 ==========
sigma = 3.0
phi_tilde = 0.019743
# ========== 小工具 ==========
def safe_filter(df: pd.DataFrame, cond_fn):
"""
对 df 应用 cond_fn(df) 的布尔筛选；若 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 get_mean_after_filter(df: pd.DataFrame, value_col: str, cond_fn):
"""
在 df 上应用 cond_fn 过滤后，取 value_col 的数值均值（忽略缺失）。
若列不存在则返回 NaN（保证稳健）。
"""
if value_col not in df.columns:
return np.nan
sub = safe_filter(df, cond_fn)
return pd.to_numeric(sub[value_col], errors="coerce").mean(skipna=True)
def gamma_from_means(cond_fn):
"""
计算：c0_bar、c1_bar、h0_bar、h1_bar 的筛选均值 → B_bar → gamma
若 B_bar <= 0 或均值缺失，返回 np.nan。
"""
c0_bar = get_mean_after_filter(c0m0full, "c0", cond_fn)
c1_bar = get_mean_after_filter(c1m1full, "c1", cond_fn)
h0_bar = get_mean_after_filter(hsf15full, "hsf15", cond_fn)
h1_bar = get_mean_after_filter(hsf18full, "hsf18", cond_fn)
# 任一均值缺失则无法计算
if any(pd.isna([c0_bar, c1_bar, h0_bar, h1_bar])):
return float("nan")
B_bar = (c0_bar ** (1 - sigma)) + (1 - sigma) * phi_tilde * (h0_bar - h1_bar)
# 对于 sigma=3，1-sigma=-2，需保证 B_bar>0 才能取幂
if not (pd.notna(B_bar) and B_bar > 0):
return float("nan")
cons1_bar = B_bar ** (1 / (1 - sigma))
gamma_val = c1_bar - cons1_bar
return float(gamma_val)
# ========== 异质性条件 ==========
# 注意：你上一条里“hsf15在5000到35000之间”应为“ic15在5000到35000之间”——这里按 ic15 实现
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, # <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
}
# ========== 批量计算并存入同名变量：gamma521 … gamma5221 ==========
_results_means = {}
for idx, cond_fn in conds.items():
name = f"gamma5{idx}1" # 末尾 1 = means-based
_results_means[name] = gamma_from_means(cond_fn)
# 可选：将结果提升为同名全局变量
globals().update(_results_means)
# 打印核对
for idx in range(2, 23):
key = f"gamma5{idx}1"
print(f"{key} = {_results_means.get(key, np.nan)}")
gamma521 = 1947.8022115045392 gamma531 = 1967.3592181966787 gamma541 = 2003.7874108596095 gamma551 = 1598.424684931106 gamma561 = 1936.6210441135872 gamma571 = 858.6784869977553 gamma581 = 2471.089904631486 gamma591 = 1453.7723798892923 gamma5101 = 1442.2229998255082 gamma5111 = 2301.3833601140914 gamma5121 = 2007.2892570366184 gamma5131 = 1618.4476362278838 gamma5141 = 2357.853554418993 gamma5151 = 1880.577649456202 gamma5161 = 1625.1534073819157 gamma5171 = 2582.726413428616 gamma5181 = 2070.0142885375976 gamma5191 = 1855.8360775174845 gamma5201 = 2768.401335332165 gamma5211 = 2298.024387545422 gamma5221 = 1737.127527122983
In [20]:
# -*- coding: utf-8 -*-
import pandas as pd
# 1) 行索引与数据
rows = [
"全样本",
"男性","女性",
"有配偶","无配偶",
"有子女","无子女",
"小于59 岁","60 岁—79 岁","80 岁及以上",
"东部","中部","西部",
"健康状况较好","健康状况中等","健康状况较差",
"较高收入","中等收入","较低收入",
"教育程度较高","教育程度中等","教育程度较低",
]
data = [
[gamma511, gamma512, gamma513],
[gamma521, gamma522, gamma523],
[gamma531, gamma532, gamma533],
[gamma541, gamma542, gamma543],
[gamma551, gamma552, gamma553],
[gamma561, gamma562, gamma563],
[gamma571, gamma572, gamma573],
[gamma581, gamma582, gamma583],
[gamma591, gamma592, gamma593],
[gamma5101, gamma5102, gamma5103],
[gamma5111, gamma5112, gamma5113],
[gamma5121, gamma5122, gamma5123],
[gamma5131, gamma5132, gamma5133],
[gamma5141, gamma5142, gamma5143],
[gamma5151, gamma5152, gamma5153],
[gamma5161, gamma5162, gamma5163],
[gamma5171, gamma5172, gamma5173],
[gamma5181, gamma5182, gamma5183],
[gamma5191, gamma5192, gamma5193],
[gamma5201, gamma5202, gamma5203],
[gamma5211, gamma5212, gamma5213],
[gamma5221, gamma5222, gamma5223],
]
# 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
# （可选）导出为 HTML 文件
# with open("表格_完全信息与最优化方法.html", "w", encoding="utf-8") as f:
# f.write(styler.to_html())
Out[20]:
| 完全信息方法 | 最优化方法 | ||
|---|---|---|---|
| 仅假设效用函数 的消费部分 | 仅假设效用函数 的健康部分 | ||
| 全样本 | 1958 | 990 | 786 |
| 男性 | 1948 | 892 | 798 |
| 女性 | 1967 | 1069 | 773 |
| 有配偶 | 2004 | 1054 | 827 |
| 无配偶 | 1598 | 592 | 549 |
| 有子女 | 1937 | 974 | 770 |
| 无子女 | 859 | 505 | 406 |
| 小于59 岁 | 2471 | 718 | 751 |
| 60 岁—79 岁 | 1454 | 1061 | 818 |
| 80 岁及以上 | 1442 | 1020 | 665 |
| 东部 | 2301 | 1167 | 726 |
| 中部 | 2007 | 1029 | 799 |
| 西部 | 1618 | 790 | 739 |
| 健康状况较好 | 2358 | 907 | 694 |
| 健康状况中等 | 1881 | 1063 | 884 |
| 健康状况较差 | 1625 | 914 | 724 |
| 较高收入 | 2583 | 945 | 793 |
| 中等收入 | 2070 | 678 | 649 |
| 较低收入 | 1856 | 1101 | 836 |
| 教育程度较高 | 2768 | 1200 | 842 |
| 教育程度中等 | 2298 | 850 | 917 |
| 教育程度较低 | 1737 | 955 | 754 |