导入数据¶
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(m0) 和 E(m1)
E_m0 = data['m0'].mean()
E_m1 = data['m1'].mean()
# 计算 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]
# 计算 gamma
data['gamma112'] = abs(data['premium2015'] - data['premium2018']) + abs(0.5 * (data['r0'] - data['r1']) * (E_m0 + E_m1)) + 0.5 * cov_c0_inv2 + 0.5 * cov_c1_inv2
gamma112= data['gamma112'].mean()
print(gamma112)
data
894.0419995170009
Out[4]:
| ID | c0 | m0 | householdID | communityID | c1 | m1 | gender | age | marriage | ... | r0 | r1 | r0adjust | r1adjust | policyintergration2015 | policyintergration2018 | district | GDPgrowthrate | urban_nbs | gamma112 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 64033321002 | 112.216 | 60.0 | 640333210 | 640333 | 123.670 | 0.0 | 0.0 | 59.0 | 1.0 | ... | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural | 704.829975 |
| 1 | 64033327002 | 1029.200 | 6000.0 | 640333270 | 640333 | 1054.764 | 21000.0 | 0.0 | 62.0 | 1.0 | ... | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural | 704.829975 |
| 2 | 64033325001 | 1062.400 | 3000.0 | 640333250 | 640333 | 78.020 | 100.0 | 1.0 | 66.0 | 1.0 | ... | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural | 704.829975 |
| 3 | 64033322001 | 592.620 | 2000.0 | 640333220 | 640333 | 859.050 | 17050.0 | 1.0 | 63.0 | 1.0 | ... | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural | 704.829975 |
| 4 | 64033330002 | 2058.400 | 2000.0 | 640333300 | 640333 | 4025.500 | 2000.0 | 1.0 | 59.0 | 1.0 | ... | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural | 704.829975 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 3882 | 89676104001 | 2315.700 | 840.0 | 896761040 | 896761 | 891.420 | 8000.0 | 1.0 | 61.0 | 1.0 | ... | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural | 753.326538 |
| 3883 | 89676114002 | 1935.145 | 300.0 | 896761140 | 896761 | 49.800 | 0.0 | 0.0 | 56.0 | 0.0 | ... | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural | 753.326538 |
| 3884 | 89676118001 | 2466.096 | 1000.0 | 896761180 | 896761 | 661.095 | 3000.0 | 0.0 | 73.0 | 0.0 | ... | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural | 753.326538 |
| 3885 | 89676115001 | 10721.940 | 800.0 | 896761150 | 896761 | 11638.260 | 2000.0 | 1.0 | 55.0 | 1.0 | ... | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural | 753.326538 |
| 3886 | 89676124001 | 268.422 | 500.0 | 896761240 | 896761 | 313.242 | 101.0 | 1.0 | 69.0 | 0.0 | ... | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural | 753.326538 |
3887 rows × 27 columns
In [6]:
#异质性男性 平衡面板
#最优化方法——消费计算
datamale = data.loc[data['gender'] == 1].copy()
import pandas as pd
import numpy as np
# 计算 E(m0) 和 E(m1)
E_m0 = datamale['m0'].mean()
E_m1 = datamale['m1'].mean()
# 计算 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]
# 计算 gamma
datamale['gamma122'] = abs(datamale['premium2015'] - datamale['premium2018']) + abs(0.5 * (datamale['r0'] - datamale['r1']) * (E_m0 + E_m1)) + 0.5 * cov_c0_inv2 + 0.5 * cov_c1_inv2
gamma122= datamale['gamma122'].mean()
print(gamma122)
795.6912224492132
In [7]:
#异质性女性 平衡面板
#最优化方法——消费计算
datafemale = data.loc[data['gender'] == 0].copy()
import pandas as pd
import numpy as np
# 计算 E(m0) 和 E(m1)
E_m0 = datafemale['m0'].mean()
E_m1 = datafemale['m1'].mean()
# 计算 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]
# 计算 gamma
datafemale['gamma132'] = abs(datafemale['premium2015'] - datafemale['premium2018']) + abs(0.5 * (datafemale['r0'] - datafemale['r1']) * (E_m0 + E_m1)) + 0.5 * cov_c0_inv2 + 0.5 * cov_c1_inv2
gamma132= datafemale['gamma132'].mean()
print(gamma132)
972.5626036980512
In [8]:
import numpy as np
import pandas as pd
# --------- 异质性条件 ---------
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 _to_num(s):
return pd.to_numeric(s, errors="coerce")
def _safe_cov(x, y):
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
# --------- 批量计算:gamma122 … gamma1222 ---------
_results = {}
for idx, cond_fn in conds.items():
# 取子样本
try:
mask = cond_fn(data)
sub = data.loc[mask].copy() if isinstance(mask, pd.Series) and len(mask)==len(data) else data.copy()
except Exception:
sub = data.copy()
if sub.empty:
_results[f"gamma1{idx}2"] = np.nan
continue
# 数值化并避免 0 的 -3 次幂产生 inf
for col in ["c0","c1","r0","r1","m0","m1","premium2015","premium2018"]:
sub[col] = _to_num(sub[col])
c0 = sub["c0"].replace(0, np.nan)
c1 = sub["c1"].replace(0, np.nan)
r0 = sub["r0"]; r1 = sub["r1"]
m0 = sub["m0"]; m1 = sub["m1"]
p15 = sub["premium2015"]; p18 = sub["premium2018"]
# 子样本的均值
E_m0 = m0.mean()
E_m1 = m1.mean()
c0_pow = c0 ** (-3)
c1_pow = c1 ** (-3)
E_c0_inv2 = c0_pow.replace([np.inf, -np.inf], np.nan).mean()
E_c1_inv2 = c1_pow.replace([np.inf, -np.inf], np.nan).mean()
# 协方差(标准化)
cov0 = cov1 = 0.0
if pd.notna(E_c0_inv2) and E_c0_inv2 != 0:
x0 = c0_pow / E_c0_inv2
y0 = (r0 - r1) * m0 + (p15 - p18)
cov0 = _safe_cov(x0, y0)
if pd.notna(E_c1_inv2) and E_c1_inv2 != 0:
x1 = c1_pow / E_c1_inv2
y1 = (r0 - r1) * m1 + (p15 - p18)
cov1 = _safe_cov(x1, y1)
# 行级 gamma,再取均值(与你示例一致)
gamma_series = (p15 - p18).abs() + (0.5 * (r0 - r1) * (E_m0 + E_m1)).abs() + 0.5*cov0 + 0.5*cov1
gamma_val = float(gamma_series.mean())
name = f"gamma1{idx}2"
_results[name] = gamma_val
globals()[name] = gamma_val # 可选:提升为同名变量
# 打印核对
for idx in range(2, 23):
key = f"gamma1{idx}2"
print(f"{key} = {_results.get(key, np.nan)}")
gamma122 = 795.6912224492132 gamma132 = 972.5626036980512 gamma142 = 968.9112657876573 gamma152 = 536.6179933332791 gamma162 = 898.4374643411936 gamma172 = 456.1738551339709 gamma182 = 638.479504002872 gamma192 = 937.0119208117399 gamma1102 = 926.4854623242564 gamma1112 = 1041.4770771710128 gamma1122 = 994.600066804255 gamma1132 = 653.0213728527 gamma1142 = 834.4207863998278 gamma1152 = 974.6068796763933 gamma1162 = 776.7776375893677 gamma1172 = 846.4128232195229 gamma1182 = 584.6022986316312 gamma1192 = 1003.8286917283358 gamma1202 = 1118.927937065343 gamma1212 = 746.2723252894853 gamma1222 = 857.9092498139793
health-based:合并数据¶
In [9]:
dataa = pd.merge(c0m0, c1m1full,on="ID", how="inner")
dataa
Out[9]:
| 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 [10]:
#计算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: 18:30:47 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[10]:
np.float64(-1.7062284276591698e-05)
In [11]:
#计算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: 18:30:48 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[11]:
np.float64(-1.3590936773055398e-05)
In [12]:
import pandas as pd
import numpy as np
# 计算 E(m0) 和 E(m1)
E_m0 = dataa['m0'].mean()
E_m1 = dataa['m1'].mean()
# 计算 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]
# 计算 gamma
dataa['gamma113'] = abs(dataa['premium2015'] - dataa['premium2018']) + abs(0.5 * (dataa['r0'] - dataa['r1']) * (E_m0 + E_m1)) + 0.5 * cov_c0_inv2 + 0.5 * cov_c1_inv2
gamma113 = dataa['gamma113'].mean()
print(gamma113)
dataa
689.4542339864197
Out[12]:
| ID | c0 | m0 | householdID | communityID | c1 | m1 | gender | age | marriage | ... | r0 | r1 | r0adjust | r1adjust | policyintergration2015 | policyintergration2018 | district | GDPgrowthrate | urban_nbs | gamma113 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 64033321002 | 112.216 | 60.0 | 640333210 | 640333 | 123.670 | 0.0 | 0.0 | 59.0 | 1.0 | ... | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural | 500.242209 |
| 1 | 64033327002 | 1029.200 | 6000.0 | 640333270 | 640333 | 1054.764 | 21000.0 | 0.0 | 62.0 | 1.0 | ... | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural | 500.242209 |
| 2 | 64033325001 | 1062.400 | 3000.0 | 640333250 | 640333 | 78.020 | 100.0 | 1.0 | 66.0 | 1.0 | ... | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural | 500.242209 |
| 3 | 64033322001 | 592.620 | 2000.0 | 640333220 | 640333 | 859.050 | 17050.0 | 1.0 | 63.0 | 1.0 | ... | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural | 500.242209 |
| 4 | 64033330002 | 2058.400 | 2000.0 | 640333300 | 640333 | 4025.500 | 2000.0 | 1.0 | 59.0 | 1.0 | ... | 0.6167 | 0.700 | 0.544192 | 0.660903 | 0.0 | 1.0 | east | 0.118005 | Rural | 500.242209 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 3882 | 89676104001 | 2315.700 | 840.0 | 896761040 | 896761 | 891.420 | 8000.0 | 1.0 | 61.0 | 1.0 | ... | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural | 548.738772 |
| 3883 | 89676114002 | 1935.145 | 300.0 | 896761140 | 896761 | 49.800 | 0.0 | 0.0 | 56.0 | 0.0 | ... | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural | 548.738772 |
| 3884 | 89676118001 | 2466.096 | 1000.0 | 896761180 | 896761 | 661.095 | 3000.0 | 0.0 | 73.0 | 0.0 | ... | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural | 548.738772 |
| 3885 | 89676115001 | 10721.940 | 800.0 | 896761150 | 896761 | 11638.260 | 2000.0 | 1.0 | 55.0 | 1.0 | ... | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural | 548.738772 |
| 3886 | 89676124001 | 268.422 | 500.0 | 896761240 | 896761 | 313.242 | 101.0 | 1.0 | 69.0 | 0.0 | ... | 0.6500 | 0.725 | 0.624462 | 0.685108 | 0.0 | 1.0 | west | 0.284050 | Rural | 548.738772 |
3887 rows × 27 columns
In [13]:
#健康based异质性探索——男性
dataamale = dataa.loc[dataa['gender'] == 1].copy()
import pandas as pd
import numpy as np
# 计算 E(m0) 和 E(m1)
E_m0 = dataamale['m0'].mean()
E_m1 = dataamale['m1'].mean()
# 计算 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]
# 计算 gamma
dataamale['gamma123'] = abs(dataamale['premium2015'] - dataamale['premium2018']) + abs(0.5 * (dataamale['r0'] - dataamale['r1']) * (E_m0 + E_m1)) + 0.5 * cov_c0_inv2 + 0.5 * cov_c1_inv2
gamma123 = dataamale['gamma123'].mean()
print(gamma123)
701.6096874269684
In [14]:
#健康based异质性探索——女性
dataafemale = dataa.loc[dataa['gender'] == 0].copy()
import pandas as pd
import numpy as np
# 计算 E(m0) 和 E(m1)
E_m0 = dataafemale['m0'].mean()
E_m1 = dataafemale['m1'].mean()
# 计算 E(c0^(-3)) 和 E(c1^(-3))
E_c0_inv2 = (dataafemale['c0']**(-3)).mean()
E_c1_inv2 = (dataafemale['c1']**(-3)).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]
# 计算 gamma
dataafemale['gamma133'] = abs(dataafemale['premium2015'] - dataafemale['premium2018']) + abs(0.5 * (dataafemale['r0'] - dataafemale['r1']) * (E_m0 + E_m1)) + 0.5 * cov_c0_inv2 + 0.5 * cov_c1_inv2
gamma133 = dataafemale['gamma133'].mean()
print(gamma133)
676.297559935399
In [15]:
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].copy()
except Exception:
pass
return df.copy()
def _to_num(s):
return pd.to_numeric(s, errors="coerce")
def _safe_cov(x, y):
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 与子样本对齐(索引或按 ID)。
支持:Series(索引=dataa.index 或 =ID)、array-like。
"""
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_explore(dataa_full, data_sub, h_m15, h_m18):
"""
你的“健康 based 探索”口径(与示例一致):
1) 子样本内均值:E_m0, E_m1;Ec0 = E[c0^(-2)], Ec1 = E[c1^(-2)]
2) 协方差:
cov0 = Cov( (PHI*h_m15)/(Ec0*r0), (r0-r1)*m0 + p15 - p18 )
cov1 = Cov( (PHI*h_m18)/(Ec1*r1), (r0-r1)*m1 + p15 - p18 )
3) 行级 gamma_i = |p15 - p18| + |0.5*(r0-r1)*(E_m0+E_m1)| + 0.5*cov0 + 0.5*cov1
返回子样本内 gamma_i 的均值
"""
# 数值化
for col in ["c0","c1","r0","r1","m0","m1","premium2015","premium2018"]:
data_sub[col] = _to_num(data_sub[col])
# 子样本均值
E_m0 = data_sub["m0"].mean()
E_m1 = data_sub["m1"].mean()
# 处理 c0/c1 的 0,避免 **(-2) 出现 inf
c0 = data_sub["c0"].replace(0, np.nan)
c1 = data_sub["c1"].replace(0, np.nan)
Ec0 = (c0 ** (-3)).replace([np.inf, -np.inf], np.nan).mean()
Ec1 = (c1 ** (-3)).replace([np.inf, -np.inf], np.nan).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)
r0 = data_sub["r0"].replace(0, np.nan)
r1 = data_sub["r1"].replace(0, np.nan)
m0 = data_sub["m0"]; m1 = data_sub["m1"]
p15 = data_sub["premium2015"]; p18 = data_sub["premium2018"]
# 协方差(Ec0/Ec1 为 0 或 NaN 时相应项置 0)
cov0 = cov1 = 0.0
if pd.notna(Ec0) and Ec0 != 0:
a0 = (PHI * h15) / (Ec0 * r0)
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)
# 行级 gamma,再取均值
gamma_series = (p15 - p18).abs() + (0.5 * (data_sub["r0"] - data_sub["r1"]) * (E_m0 + E_m1)).abs() + 0.5*cov0 + 0.5*cov1
return float(gamma_series.mean())
# ========== 异质性条件 ==========
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
}
# ========== 批量计算:gamma123 … gamma1223 ==========
_results = {}
for idx, cond_fn in conds.items():
sub = safe_filter(dataa, cond_fn) # 按你的示例,health-based 用 dataa 作为主样本
if sub.empty:
val = np.nan
else:
val = compute_gamma_health_explore(dataa_full=dataa, data_sub=sub, h_m15=h_m15, h_m18=h_m18)
name = f"gamma1{idx}3"
_results[name] = val
globals()[name] = val # 可选:注册为同名变量
# 打印核对
for idx in range(2, 23):
key = f"gamma1{idx}3"
print(f"{key} = {_results.get(key, np.nan)}")
gamma123 = 701.6096874269684 gamma133 = 676.297559935399 gamma143 = 742.1281289436823 gamma153 = 493.6350274056347 gamma163 = 694.3024694656491 gamma173 = 356.45453692191603 gamma183 = 671.0503320289147 gamma193 = 693.9424656194997 gamma1103 = 571.7774025800886 gamma1113 = 600.9391876768528 gamma1123 = 763.7401963175519 gamma1133 = 601.8149627531205 gamma1143 = 621.8740168169526 gamma1153 = 795.3126953827027 gamma1163 = 586.7692087079621 gamma1173 = 694.1307975452706 gamma1183 = 555.6822076078612 gamma1193 = 739.0296290542765 gamma1203 = 760.9427888447271 gamma1213 = 813.6137702363884 gamma1223 = 656.7930492144557
用完全信息法求解¶
In [16]:
d1 = pd.merge(c0m0, c1m1, on="ID", how="inner")
d2 = pd.merge(d1, hsf15, on="ID", how="inner")
d3 = pd.merge(d2, hsf18, on="ID", how="inner")
d3
Out[16]:
| ID | c0 | m0 | c1 | m1 | hsf15 | hsf18 | |
|---|---|---|---|---|---|---|---|
| 0 | 64033321002 | 112.216 | 60.0 | 123.670 | 0.0 | 54.436016 | 10.551840 |
| 1 | 64033327002 | 1029.200 | 6000.0 | 1054.764 | 21000.0 | 28.096293 | 39.754467 |
| 2 | 64033322001 | 592.620 | 2000.0 | 859.050 | 17050.0 | 21.303569 | 59.803845 |
| 3 | 64033330002 | 2058.400 | 2000.0 | 4025.500 | 2000.0 | 34.523055 | 68.626626 |
| 4 | 64033341001 | 5436.500 | 500.0 | 1806.578 | 2500.0 | 50.384267 | 78.657775 |
| ... | ... | ... | ... | ... | ... | ... | ... |
| 3750 | 89676104001 | 2315.700 | 840.0 | 891.420 | 8000.0 | 27.522919 | 42.058152 |
| 3751 | 89676114002 | 1935.145 | 300.0 | 49.800 | 0.0 | 35.235805 | 15.917436 |
| 3752 | 89676118001 | 2466.096 | 1000.0 | 661.095 | 3000.0 | 32.251424 | 54.653869 |
| 3753 | 89676115001 | 10721.940 | 800.0 | 11638.260 | 2000.0 | 32.757347 | 73.821948 |
| 3754 | 89676124001 | 268.422 | 500.0 | 313.242 | 101.0 | 34.706036 | 95.825352 |
3755 rows × 7 columns
In [17]:
import numpy as np
import pandas as pd
# 参数
sigma = 3.0
phi_tilde = 0.019743
d3["B_bar"] = (d3["c0"]**(1 - sigma)) + (1 - sigma) * phi_tilde * (d3["hsf15"] - d3["hsf18"])
d3["gamma111"] = d3["c1"] - d3["B_bar"]**(1 / (1 - sigma))
gamma111= d3["gamma111"].mean()
print(gamma111)
d3
1996.3867642406108
Out[17]:
| ID | c0 | m0 | c1 | m1 | hsf15 | hsf18 | B_bar | gamma111 | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 64033321002 | 112.216 | 60.0 | 123.670 | 0.0 | 54.436016 | 10.551840 | -1.732731 | NaN |
| 1 | 64033327002 | 1029.200 | 6000.0 | 1054.764 | 21000.0 | 28.096293 | 39.754467 | 0.460336 | 1053.290118 |
| 2 | 64033322001 | 592.620 | 2000.0 | 859.050 | 17050.0 | 21.303569 | 59.803845 | 1.520225 | 858.238953 |
| 3 | 64033330002 | 2058.400 | 2000.0 | 4025.500 | 2000.0 | 34.523055 | 68.626626 | 1.346614 | 4024.638256 |
| 4 | 64033341001 | 5436.500 | 500.0 | 1806.578 | 2500.0 | 50.384267 | 78.657775 | 1.116408 | 1805.631570 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 3750 | 89676104001 | 2315.700 | 840.0 | 891.420 | 8000.0 | 27.522919 | 42.058152 | 0.573938 | 890.100020 |
| 3751 | 89676114002 | 1935.145 | 300.0 | 49.800 | 0.0 | 35.235805 | 15.917436 | -0.762805 | NaN |
| 3752 | 89676118001 | 2466.096 | 1000.0 | 661.095 | 3000.0 | 32.251424 | 54.653869 | 0.884583 | 660.031762 |
| 3753 | 89676115001 | 10721.940 | 800.0 | 11638.260 | 2000.0 | 32.757347 | 73.821948 | 1.621477 | 11637.474684 |
| 3754 | 89676124001 | 268.422 | 500.0 | 313.242 | 101.0 | 34.706036 | 95.825352 | 2.413371 | 312.598293 |
3755 rows × 9 columns
In [18]:
#异质性男性
import numpy as np
import pandas as pd
d4 = pd.merge(d1, hsf18full, on="ID", how="inner")
dmale = d4.loc[d4['gender'] == 1].copy()
# 参数
sigma = 3.0
phi_tilde = 0.019743
dmale["B_bar"] = (dmale["c0"]**(1 - sigma)) + (1 - sigma) * phi_tilde * (dmale["hsf15"] - dmale["hsf18"])
dmale["gamma121"] = dmale["c1"] - dmale["B_bar"]**(1 / (1 - sigma))
gamma121= dmale["gamma121"].mean()
print(gamma121)
1968.4893177925992
In [19]:
#异质性女性
import numpy as np
import pandas as pd
d4 = pd.merge(d1, hsf18full, on="ID", how="inner")
dfemale = d4.loc[d4['gender'] == 0].copy()
# 参数
sigma = 3.0
phi_tilde = 0.019743
dfemale["B_bar"] = (dfemale["c0"]**(1 - sigma)) + (1 - sigma) * phi_tilde * (dfemale["hsf15"] - dfemale["hsf18"])
dfemale["gamma131"] = dfemale["c1"] - dfemale["B_bar"]**(1 / (1 - sigma))
gamma131= dfemale["gamma131"].mean()
print(gamma131)
2030.5047297378765
In [20]:
import numpy as np
import pandas as pd
# ===== 1) 构造 d4 =====
# 这里假定 d1 已经包含了 hsf15、c0、c1、gender 等列
d4 = pd.merge(d1, hsf18full, on="ID", how="inner")
# 若某些分组变量是字符串(比如 district),确保为字符串
if "district" in d4.columns:
d4["district"] = d4["district"].astype(str)
# ===== 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, # 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"].str.lower().eq("east"), # east
12: lambda d: d["district"].str.lower().eq("middle"), # middle
13: lambda d: d["district"].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) 核心计算(逐行算,再对子样本取均值) =====
sigma = 3.0
phi_tilde = 0.019743
power = 1.0 - sigma # = -2
inv_power = 1.0 / power # = -0.5
def compute_gamma_subset_like_yours(df_sub: pd.DataFrame) -> float:
if df_sub.empty:
return float("nan")
# 不做 base>0 保护
B_bar = (df_sub["c0"] ** power) + power * phi_tilde * (df_sub["hsf15"] - df_sub["hsf18"])
gamma_i = df_sub["c1"] - (B_bar ** inv_power)
return float(gamma_i.mean()) # mean 默认会跳过 NaN
# ===== 4) 批量生成 gamma121 … gamma1221 =====
results = {}
for idx, cond_fn in conds.items():
# 如果条件里用到的列不存在,会抛错;这里捕获后把该组记为 NaN(避免改变其它组结果)
try:
mask = cond_fn(d4)
sub = d4.loc[mask].copy() if isinstance(mask, pd.Series) and len(mask) == len(d4) else d4.copy()
except Exception:
sub = pd.DataFrame(columns=d4.columns) # 列缺失时该组置空
name = f"gamma1{idx}1"
results[name] = compute_gamma_subset_like_yours(sub)
globals()[name] = results[name] # 可选:注册为同名变量
# 打印核对
for idx in range(2, 23):
key = f"gamma1{idx}1"
print(f"{key} = {results.get(key, np.nan)}")
gamma121 = 1968.4893177925992 gamma131 = 2030.5047297378765 gamma141 = 2091.756836228619 gamma151 = 1516.9958081205616 gamma161 = 2009.6788446394628 gamma171 = 752.4025142905097 gamma181 = 2513.485763337119 gamma191 = 1492.2099449302036 gamma1101 = 963.7132207091026 gamma1111 = 2318.3816169376114 gamma1121 = 2056.9974605815864 gamma1131 = 1614.5882274330595 gamma1141 = 2312.716443873066 gamma1151 = 1938.9117698957555 gamma1161 = 1559.7312849130187 gamma1171 = 2566.5778502883627 gamma1181 = 2138.1528481056553 gamma1191 = 1882.5176488438392 gamma1201 = 2677.6344623042633 gamma1211 = 2276.608300019149 gamma1221 = 1740.557558427637
In [21]:
# -*- coding: utf-8 -*-
import pandas as pd
# 1) 行索引与数据
rows = [
"全样本",
"男性","女性",
"有配偶","无配偶",
"有子女","无子女",
"小于59 岁","60 岁—79 岁","80 岁及以上",
"东部","中部","西部",
"健康状况较好","健康状况中等","健康状况较差",
"较高收入","中等收入","较低收入",
"教育程度较高","教育程度中等","教育程度较低",
]
data = [
[gamma111, gamma112, gamma113],
[gamma121, gamma122, gamma123],
[gamma131, gamma132, gamma133],
[gamma141, gamma142, gamma143],
[gamma151, gamma152, gamma153],
[gamma161, gamma162, gamma163],
[gamma171, gamma172, gamma173],
[gamma181, gamma182, gamma183],
[gamma191, gamma192, gamma193],
[gamma1101, gamma1102, gamma1103],
[gamma1111, gamma1112, gamma1113],
[gamma1121, gamma1122, gamma1123],
[gamma1131, gamma1132, gamma1133],
[gamma1141, gamma1142, gamma1143],
[gamma1151, gamma1152, gamma1153],
[gamma1161, gamma1162, gamma1163],
[gamma1171, gamma1172, gamma1173],
[gamma1181, gamma1182, gamma1183],
[gamma1191, gamma1192, gamma1193],
[gamma1201, gamma1202, gamma1203],
[gamma1211, gamma1212, gamma1213],
[gamma1221, gamma1222, gamma1223],
]
# 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[21]:
| 完全信息方法 | 最优化方法 | ||
|---|---|---|---|
| 仅假设效用函数 的消费部分 | 仅假设效用函数 的健康部分 | ||
| 全样本 | 1996 | 894 | 689 |
| 男性 | 1968 | 796 | 702 |
| 女性 | 2031 | 973 | 676 |
| 有配偶 | 2092 | 969 | 742 |
| 无配偶 | 1517 | 537 | 494 |
| 有子女 | 2010 | 898 | 694 |
| 无子女 | 752 | 456 | 356 |
| 小于59 岁 | 2513 | 638 | 671 |
| 60 岁—79 岁 | 1492 | 937 | 694 |
| 80 岁及以上 | 964 | 926 | 572 |
| 东部 | 2318 | 1041 | 601 |
| 中部 | 2057 | 995 | 764 |
| 西部 | 1615 | 653 | 602 |
| 健康状况较好 | 2313 | 834 | 622 |
| 健康状况中等 | 1939 | 975 | 795 |
| 健康状况较差 | 1560 | 777 | 587 |
| 较高收入 | 2567 | 846 | 694 |
| 中等收入 | 2138 | 585 | 556 |
| 较低收入 | 1883 | 1004 | 739 |
| 教育程度较高 | 2678 | 1119 | 761 |
| 教育程度中等 | 2277 | 746 | 814 |
| 教育程度较低 | 1741 | 858 | 657 |