导入数据

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')  

#删除城市样本
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)]

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']

#删去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：合并数据

data = pd.merge(c0m0, c1m1full,on="ID", how="inner")
data

	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

#最优化方法——消费计算
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['gamma412'] = 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
gamma412= data['gamma412'].mean()
print(gamma412)
data

868.3633602368919

	ID	c0	m0	householdID	communityID	c1	m1	gender	age	marriage	...	r0	r1	r0adjust	r1adjust	policyintergration2015	policyintergration2018	district	GDPgrowthrate	urban_nbs	gamma412
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	673.268989
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	673.268989
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	673.268989
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	673.268989
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	673.268989
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
3258	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	721.976891
3259	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	721.976891
3260	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	721.976891
3261	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	721.976891
3262	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	721.976891

3263 rows × 27 columns

#异质性男性 平衡面板
#最优化方法——消费计算 
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['gamma422'] = 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
gamma422= datamale['gamma422'].mean()
print(gamma422)

694.7652383900032

#异质性女性 平衡面板
#最优化方法——消费计算 
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['gamma432'] = 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
gamma432= datafemale['gamma432'].mean()
print(gamma432)

977.057580096069

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 的 -2 次幂产生 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"gamma4{idx}2"
    _results[name] = gamma_val
    globals()[name] = gamma_val  # 可选：提升为同名变量

# 打印核对
for idx in range(2, 23):
    key = f"gamma4{idx}2"
    print(f"{key} = {_results.get(key, np.nan)}")

gamma422 = 694.7652383900032
gamma432 = 977.057580096069
gamma442 = 928.9212638897537
gamma452 = 524.1607537999155
gamma462 = 870.8773550140032
gamma472 = 503.0717633593451
gamma482 = 622.4214532204546
gamma492 = 873.8549825967619
gamma4102 = 912.6547535792567
gamma4112 = 908.1031871941158
gamma4122 = 935.7939737706656
gamma4132 = 725.8448703364339
gamma4142 = 787.6652592920365
gamma4152 = 922.5700216388028
gamma4162 = 804.7370749835521
gamma4172 = 740.3300133547918
gamma4182 = 538.9672255756258
gamma4192 = 990.5126460339013
gamma4202 = 1147.9738699490313
gamma4212 = 689.5772417336428
gamma4222 = 818.7179387434956

health-based：合并数据

dataa = pd.merge(c0m0, c1m1full,on="ID", how="inner")
dataa

	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

#计算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:                        21:55:28   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.

np.float64(-1.72380438213464e-05)

#计算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:                        21:55:28   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.

np.float64(-1.3328744662871004e-05)

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['gamma413'] = 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
gamma413 = dataa['gamma413'].mean()
print(gamma413)
dataa

696.4065068044124

	ID	c0	m0	householdID	communityID	c1	m1	gender	age	marriage	...	r0	r1	r0adjust	r1adjust	policyintergration2015	policyintergration2018	district	GDPgrowthrate	urban_nbs	gamma413
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	501.312136
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	501.312136
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	501.312136
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	501.312136
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	501.312136
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
3258	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	550.020038
3259	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	550.020038
3260	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	550.020038
3261	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	550.020038
3262	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	550.020038

3263 rows × 27 columns

#健康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['gamma423'] = 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
gamma423 = dataamale['gamma423'].mean()
print(gamma423)

713.2089647232881

#健康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['gamma433'] = 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
gamma433 = dataafemale['gamma433'].mean()
print(gamma433)

686.9411944892615

import numpy as np
import pandas as pd

PHI = 0.019743

# ========== 工具函数 ==========
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):
    """
      1) 子样本内均值：E_m0, E_m1；Ec0 = E[c0^(-3)], Ec1 = E[c1^(-3)]
      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
    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) 
    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"gamma4{idx}3"
    _results[name] = val
    globals()[name] = val 

# 打印核对
for idx in range(2, 23):
    key = f"gamma4{idx}3"
    print(f"{key} = {_results.get(key, np.nan)}")

gamma423 = 713.2089647232881
gamma433 = 686.9411944892615
gamma443 = 755.2935186124944
gamma453 = 495.3709187196489
gamma463 = 701.0600646792204
gamma473 = 384.02320419050864
gamma483 = 692.6520656288419
gamma493 = 680.4721585389249
gamma4103 = 554.6613086282168
gamma4113 = 538.9687740426275
gamma4123 = 790.6987379683895
gamma4133 = 637.7637835990973
gamma4143 = 623.7171886986312
gamma4153 = 812.5947344792033
gamma4163 = 605.8677243993955
gamma4173 = 634.2296693498457
gamma4183 = 555.12008641023
gamma4193 = 754.2440697210108
gamma4203 = 779.6733659370648
gamma4213 = 835.7746510292079
gamma4223 = 657.4253712423545

用完全信息法求解

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

	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
...	...	...	...	...	...	...	...
3145	89676104001	2315.700	840.0	891.420	8000.0	27.522919	42.058152
3146	89676114002	1935.145	300.0	49.800	0.0	35.235805	15.917436
3147	89676118001	2466.096	1000.0	661.095	3000.0	32.251424	54.653869
3148	89676115001	10721.940	800.0	11638.260	2000.0	32.757347	73.821948
3149	89676124001	268.422	500.0	313.242	101.0	34.706036	95.825352

3150 rows × 7 columns

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["gamma411"] = d3["c1"] - d3["B_bar"]**(1 / (1 - sigma))

gamma411= d3["gamma411"].mean()
print(gamma411)
d3

1966.4793432831154

	ID	c0	m0	c1	m1	hsf15	hsf18	B_bar	gamma411
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
...	...	...	...	...	...	...	...	...	...
3145	89676104001	2315.700	840.0	891.420	8000.0	27.522919	42.058152	0.573938	890.100020
3146	89676114002	1935.145	300.0	49.800	0.0	35.235805	15.917436	-0.762805	NaN
3147	89676118001	2466.096	1000.0	661.095	3000.0	32.251424	54.653869	0.884583	660.031762
3148	89676115001	10721.940	800.0	11638.260	2000.0	32.757347	73.821948	1.621477	11637.474684
3149	89676124001	268.422	500.0	313.242	101.0	34.706036	95.825352	2.413371	312.598293

3150 rows × 9 columns

#异质性男性
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["gamma421"] = dmale["c1"] - dmale["B_bar"]**(1 / (1 - sigma))

gamma421= dmale["gamma421"].mean()
print(gamma421)

1926.768109895933

#异质性女性
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["gamma431"] = dfemale["c1"] - dfemale["B_bar"]**(1 / (1 - sigma))

gamma431= dfemale["gamma431"].mean()
print(gamma431)

2012.4698507593214

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"gamma4{idx}1"
    results[name] = compute_gamma_subset_like_yours(sub)
    globals()[name] = results[name] 

# 打印核对
for idx in range(2, 23):
    key = f"gamma4{idx}1"
    print(f"{key} = {results.get(key, np.nan)}")

gamma421 = 1926.768109895933
gamma431 = 2012.4698507593214
gamma441 = 2066.4616597642485
gamma451 = 1454.7471914412529
gamma461 = 1978.5411036233604
gamma471 = 783.4967750642056
gamma481 = 2498.9166338774003
gamma491 = 1449.6459516847665
gamma4101 = 1052.1186704357992
gamma4111 = 2262.993445095227
gamma4121 = 2027.3414098113271
gamma4131 = 1580.7204444612387
gamma4141 = 2301.6770212805272
gamma4151 = 1909.327378044254
gamma4161 = 1481.321882864932
gamma4171 = 2480.90139971734
gamma4181 = 2063.6424961812704
gamma4191 = 1873.807932658593
gamma4201 = 2643.282946856107
gamma4211 = 2279.502473938127
gamma4221 = 1705.0261814509524

# -*- coding: utf-8 -*-
import pandas as pd

# 1) 行索引与数据
rows = [
    "全样本",
    "男性","女性",
    "有配偶","无配偶",
    "有子女","无子女",
    "小于59 岁","60 岁—79 岁","80 岁及以上",
    "东部","中部","西部",
    "健康状况较好","健康状况中等","健康状况较差",
    "较高收入","中等收入","较低收入",
    "教育程度较高","教育程度中等","教育程度较低",
]

data = [
[gamma411, gamma412, gamma413],
[gamma421, gamma422, gamma423],
[gamma431, gamma432, gamma433],
[gamma441, gamma442, gamma443],
[gamma451, gamma452, gamma453],
[gamma461, gamma462, gamma463],
[gamma471, gamma472, gamma473],
[gamma481, gamma482, gamma483],
[gamma491, gamma492, gamma493],
[gamma4101, gamma4102, gamma4103],
[gamma4111, gamma4112, gamma4113],
[gamma4121, gamma4122, gamma4123],
[gamma4131, gamma4132, gamma4133],
[gamma4141, gamma4142, gamma4143],
[gamma4151, gamma4152, gamma4153],
[gamma4161, gamma4162, gamma4163],
[gamma4171, gamma4172, gamma4173],
[gamma4181, gamma4182, gamma4183],
[gamma4191, gamma4192, gamma4193],
[gamma4201, gamma4202, gamma4203],
[gamma4211, gamma4212, gamma4213],
[gamma4221, gamma4222, gamma4223],
]

# 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

	完全信息方法	最优化方法
		仅假设效用函数 的消费部分	仅假设效用函数 的健康部分
全样本	1966	868	696
男性	1927	695	713
女性	2012	977	687
有配偶	2066	929	755
无配偶	1455	524	495
有子女	1979	871	701
无子女	783	503	384
小于59 岁	2499	622	693
60 岁—79 岁	1450	874	680
80 岁及以上	1052	913	555
东部	2263	908	539
中部	2027	936	791
西部	1581	726	638
健康状况较好	2302	788	624
健康状况中等	1909	923	813
健康状况较差	1481	805	606
较高收入	2481	740	634
中等收入	2064	539	555
较低收入	1874	991	754
教育程度较高	2643	1148	780
教育程度中等	2280	690	836
教育程度较低	1705	819	657