导入数据¶
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)]
#删除城市样本
c0m0full = c0m0full[c0m0full['urban_nbs'] != 'Urban']
#确保15年未整合，18年整合了
c0m0full = c0m0full[(c0m0full['policyintergration2015']==0.0) & (c0m0full['policyintergration2018']==1.0)]
#删除城市样本
hsf15full = hsf15full[hsf15full['urban_nbs'] != 'Urban']
#确保15年未整合，18年整合了
hsf15full = hsf15full[(hsf15full['policyintergration2015']==0.0) & (hsf15full['policyintergration2018']==1.0)]
#删除城市样本
hsf18full = hsf18full[hsf18full['urban_nbs'] != 'Urban']
#确保15年未整合，18年整合了
hsf18full = hsf18full[(hsf18full['policyintergration2015']==0.0) & (hsf18full['policyintergration2018']==1.0)]
In [3]:
import numpy as np
import pandas as pd
# 1) 只保留需要列，去除缺失，按 city 去重（同城增速一致）
dfu = (city_gender_age_premium_ratio_district15[['city','GDPgrowthrate']]
.dropna(subset=['GDPgrowthrate'])
.drop_duplicates(subset=['city']))
# 2) 按增速从高到低排序
df_sorted = dfu.sort_values('GDPgrowthrate', ascending=False).reset_index(drop=True)
# 3) 90%分位（排名法）：正数10%位置的阈值 = 从上往下第 ceil(0.1*N) 个
n = len(df_sorted)
k = max(1, int(np.ceil(0.10 * n)))
threshold = df_sorted.loc[k - 1, 'GDPgrowthrate']
# 4) 列出严格高于阈值的城市
top_cities = df_sorted.loc[df_sorted['GDPgrowthrate'] > threshold, 'city'].tolist()
print("90%分位（排名法，正数10%位置）k =", k, "阈值 GDPgrowthrate =", threshold)
print("高于该阈值的城市：", top_cities)
90%分位（排名法，正数10%位置）k = 13 阈值 GDPgrowthrate = 0.384969765991486 高于该阈值的城市： ['Yulin', 'Qiannan', 'Xinzhou', 'Chuxiong', 'Zhangzhou', 'Ganzhou', 'Chengdu', 'Jiujiang', 'Fuzhou', 'Chaohu', 'Fuyang', 'Zhengzhou']
In [4]:
#删去GDP增速过快的城市
c0m0full = c0m0full[c0m0full['GDPgrowthrate']<=threshold]
c0m0= c0m0full[['ID', 'c0','m0']]
c1m1full = c1m1full[c1m1full['GDPgrowthrate']<=threshold]
c1m1= c1m1full[['ID', 'c1','m1']]
hsf15full = hsf15full[hsf15full['GDPgrowthrate']<=threshold]
hsf15= hsf15full[['ID', 'hsf15']]
hsf18full = hsf18full[hsf18full['GDPgrowthrate']<=threshold]
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 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 3258 | 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 |
| 3259 | 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 |
| 3260 | 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 |
| 3261 | 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 |
| 3262 | 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 |
3263 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混合截面）
gamma812=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(gamma812)
Out[6]:
954.7280836011141
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]
gamma822=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(gamma822)
Out[7]:
786.9805456609977
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]
gamma832=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(gamma832)
Out[8]:
1058.060165337711
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
}
# —— 批量计算，每个条件生成一个 gamma7{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'gamma8{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'gamma8{idx}2'
print(f'{key} = {_results.get(key, np.nan)}')
gamma822 = 786.9805456609977 gamma832 = 1058.060165337711 gamma842 = 1003.8021842982591 gamma852 = 583.985827155903 gamma862 = 939.9856923166444 gamma872 = 530.9929611436345 gamma882 = 685.4535815104775 gamma892 = 994.8066281751774 gamma8102 = 1010.1010661694434 gamma8112 = 1039.7627672008864 gamma8122 = 961.3615452245034 gamma8132 = 831.8667377013901 gamma8142 = 861.0847272264871 gamma8152 = 991.3379058580065 gamma8162 = 929.231391748959 gamma8172 = 872.7998986960389 gamma8182 = 618.7220358460878 gamma8192 = 1075.3279431279848 gamma8202 = 1236.8656898156946 gamma8212 = 792.5012580950643 gamma8222 = 901.4947694070074
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 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 3258 | 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 |
| 3259 | 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 |
| 3260 | 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 |
| 3261 | 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 |
| 3262 | 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 |
3263 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.124
Date: Mon, 29 Dec 2025 Prob (F-statistic): 0.289
Time: 18:14:05 Log-Likelihood: -10028.
No. Observations: 2607 AIC: 2.006e+04
Df Residuals: 2605 BIC: 2.007e+04
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 33.6034 0.237 141.906 0.000 33.139 34.068
m0 -1.724e-05 1.63e-05 -1.060 0.289 -4.91e-05 1.46e-05
==============================================================================
Omnibus: 37.273 Durbin-Watson: 1.858
Prob(Omnibus): 0.000 Jarque-Bera (JB): 28.788
Skew: -0.170 Prob(JB): 5.61e-07
Kurtosis: 2.613 Cond. No. 1.55e+04
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.55e+04. This might indicate that there are
strong multicollinearity or other numerical problems.
Out[11]:
np.float64(-1.72380438213464e-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.259
Date: Mon, 29 Dec 2025 Prob (F-statistic): 0.262
Time: 18:14:06 Log-Likelihood: -22129.
No. Observations: 4727 AIC: 4.426e+04
Df Residuals: 4725 BIC: 4.427e+04
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 52.1728 0.391 133.428 0.000 51.406 52.939
m1 -1.333e-05 1.19e-05 -1.122 0.262 -3.66e-05 9.95e-06
==============================================================================
Omnibus: 7734.392 Durbin-Watson: 1.732
Prob(Omnibus): 0.000 Jarque-Bera (JB): 372.216
Skew: -0.275 Prob(JB): 1.49e-81
Kurtosis: 1.740 Cond. No. 3.39e+04
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 3.39e+04. This might indicate that there are
strong multicollinearity or other numerical problems.
Out[12]:
np.float64(-1.3328744662871004e-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]
gamma813=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(gamma813)
Out[13]:
782.7712301686346
In [14]:
#健康based异质性探索——男性
dataamale=dataa[dataa['gender'] == 1]
import pandas as pd
import numpy as np
# 计算 E(c0^(-2)) 和 E(c1^(-2))
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混合截面）
gamma823=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(gamma823)
Out[14]:
805.4242719942827
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) 条件（编号 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, # <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) 批量计算：gamma423 … gamma4223 ==================
_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"gamma8{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"gamma8{idx}3"
print(f"{key} = {_results_h.get(key, np.nan)}")
gamma823 = 805.4242719942827 gamma833 = 767.9437797309033 gamma843 = 830.1744390209999 gamma853 = 555.1959920756364 gamma863 = 770.1684019818616 gamma873 = 411.94440197479804 gamma883 = 755.6841939188649 gamma893 = 801.4238041173403 gamma8103 = 652.1076212184036 gamma8113 = 670.6283540493981 gamma8123 = 816.2663094222273 gamma8133 = 743.7856509640535 gamma8143 = 697.1366566330818 gamma8153 = 881.3626186984069 gamma8163 = 730.3620411648025 gamma8173 = 766.6995546910928 gamma8183 = 634.8748966806921 gamma8193 = 839.0593668150944 gamma8203 = 868.5651858037281 gamma8213 = 938.6986673906292 gamma8223 = 740.2022019058663
用完全信息法求解¶
In [16]:
import numpy as np
import pandas as pd
# 参数
sigmamale = 3.0
phi_tilde = 0.019743
# 取各列的均值（忽略缺失）
c0_bar = pd.to_numeric(c0m0full["c0"], errors="coerce").mean(skipna=True)
c1_bar = pd.to_numeric(c1m1full["c1"], errors="coerce").mean(skipna=True)
h0_bar = pd.to_numeric(hsf15full["hsf15"], errors="coerce").mean(skipna=True)
h1_bar = pd.to_numeric(hsf18full["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))
gamma811 = c1_bar - cons1_bar
print(gamma811)
1927.543007125511
In [17]:
#异质性计算男性
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))
gamma821 = c1_bar - cons1_bar
print(gamma821)
1922.2866178040815
In [18]:
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)
# ========== 异质性条件（编号 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, # <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
}
# ========== 批量计算并存入同名变量：gamma421 … gamma4221 ==========
_results_means = {}
for idx, cond_fn in conds.items():
name = f"gamma8{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"gamma8{idx}1"
print(f"{key} = {_results_means.get(key, np.nan)}")
gamma821 = 1922.2866178040815 gamma831 = 1933.0937413047204 gamma841 = 1977.03035451382 gamma851 = 1535.7053945281282 gamma861 = 1901.303458439417 gamma871 = 938.354330523831 gamma881 = 2438.4528777720298 gamma891 = 1433.4701563394885 gamma8101 = 1404.6022327483054 gamma8111 = 2271.950229493393 gamma8121 = 1974.139732933464 gamma8131 = 1582.7674094223478 gamma8141 = 2308.73491240512 gamma8151 = 1860.6803872618227 gamma8161 = 1584.6586640366686 gamma8171 = 2404.580782022482 gamma8181 = 2023.0556849982147 gamma8191 = 1843.9424748220827 gamma8201 = 2679.300448644135 gamma8211 = 2312.4095212690004 gamma8221 = 1693.2846470465943
In [19]:
# -*- coding: utf-8 -*-
import pandas as pd
# 1) 行索引与数据
rows = [
"全样本",
"男性","女性",
"有配偶","无配偶",
"有子女","无子女",
"小于59 岁","60 岁—79 岁","80 岁及以上",
"东部","中部","西部",
"健康状况较好","健康状况中等","健康状况较差",
"较高收入","中等收入","较低收入",
"教育程度较高","教育程度中等","教育程度较低",
]
data = [
[gamma811, gamma812, gamma813],
[gamma821, gamma822, gamma823],
[gamma831, gamma832, gamma833],
[gamma841, gamma842, gamma843],
[gamma851, gamma852, gamma853],
[gamma861, gamma862, gamma863],
[gamma871, gamma872, gamma873],
[gamma881, gamma882, gamma883],
[gamma891, gamma892, gamma893],
[gamma8101, gamma8102, gamma8103],
[gamma8111, gamma8112, gamma8113],
[gamma8121, gamma8122, gamma8123],
[gamma8131, gamma8132, gamma8133],
[gamma8141, gamma8142, gamma8143],
[gamma8151, gamma8152, gamma8153],
[gamma8161, gamma8162, gamma8163],
[gamma8171, gamma8172, gamma8173],
[gamma8181, gamma8182, gamma8183],
[gamma8191, gamma8192, gamma8193],
[gamma8201, gamma8202, gamma8203],
[gamma8211, gamma8212, gamma8213],
[gamma8221, gamma8222, gamma8223],
]
# 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[19]:
| 完全信息方法 | 最优化方法 | ||
|---|---|---|---|
| 仅假设效用函数 的消费部分 | 仅假设效用函数 的健康部分 | ||
| 全样本 | 1928 | 955 | 783 |
| 男性 | 1922 | 787 | 805 |
| 女性 | 1933 | 1058 | 768 |
| 有配偶 | 1977 | 1004 | 830 |
| 无配偶 | 1536 | 584 | 555 |
| 有子女 | 1901 | 940 | 770 |
| 无子女 | 938 | 531 | 412 |
| 小于59 岁 | 2438 | 685 | 756 |
| 60 岁—79 岁 | 1433 | 995 | 801 |
| 80 岁及以上 | 1405 | 1010 | 652 |
| 东部 | 2272 | 1040 | 671 |
| 中部 | 1974 | 961 | 816 |
| 西部 | 1583 | 832 | 744 |
| 健康状况较好 | 2309 | 861 | 697 |
| 健康状况中等 | 1861 | 991 | 881 |
| 健康状况较差 | 1585 | 929 | 730 |
| 较高收入 | 2405 | 873 | 767 |
| 中等收入 | 2023 | 619 | 635 |
| 较低收入 | 1844 | 1075 | 839 |
| 教育程度较高 | 2679 | 1237 | 869 |
| 教育程度中等 | 2312 | 793 | 939 |
| 教育程度较低 | 1693 | 901 | 740 |