导入数据

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)]
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：合并数据

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
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
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

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

894.0419995170009

	ID	c0	m0	householdID	communityID	c1	m1	gender	age	marriage	...	r0	r1	r0adjust	r1adjust	policyintergration2015	policyintergration2018	district	GDPgrowthrate	urban_nbs	gamma212
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

#异质性男性 平衡面板
#最优化方法——消费计算 
datamale = data.loc[data['gender'] == 1].copy()
sigamamale = 2.6
import pandas as pd
import numpy as np

# 计算 E(m0) 和 E(m1)
E_m0 = datamale['m0'].mean()
E_m1 = datamale['m1'].mean()

# 计算 E(c0^(-sigma)) 和 E(c1^(-sigma))
E_c0_inv2 = (datamale['c0']**(-sigamamale)).mean()
E_c1_inv2 = (datamale['c1']**(-sigamamale)).mean()

# 计算协方差
cov_c0_inv2 = np.cov(datamale['c0']**(-sigamamale) / E_c0_inv2, (datamale['r0'] - datamale['r1']) * datamale['m0'] + datamale['premium2015'] - datamale['premium2018'])[0, 1]
cov_c1_inv2 = np.cov(datamale['c1']**(-sigamamale) / E_c1_inv2, (datamale['r0'] - datamale['r1']) * datamale['m1'] + datamale['premium2015'] - datamale['premium2018'])[0, 1]

# 计算 gamma
datamale['gamma222'] = 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
gamma222= datamale['gamma222'].mean()
print(gamma222)

806.2642603893461

#异质性女性 平衡面板
#最优化方法——消费计算 
datafemale = data.loc[data['gender'] == 0].copy()
sigamafemale = 3.4
import pandas as pd
import numpy as np

# 计算 E(m0) 和 E(m1)
E_m0 = datafemale['m0'].mean()
E_m1 = datafemale['m1'].mean()

# 计算 E(c0^(-sigma)) 和 E(c1^(-sigma))
E_c0_inv2 = (datafemale['c0']**(-sigamafemale)).mean()
E_c1_inv2 = (datafemale['c1']**(-sigamafemale)).mean()

# 计算协方差
cov_c0_inv2 = np.cov(datafemale['c0']**(-sigamafemale) / E_c0_inv2, (datafemale['r0'] - datafemale['r1']) * datafemale['m0'] + datafemale['premium2015'] - datafemale['premium2018'])[0, 1]
cov_c1_inv2 = np.cov(datafemale['c1']**(-sigamafemale) / E_c1_inv2, (datafemale['r0'] - datafemale['r1']) * datafemale['m1'] + datafemale['premium2015'] - datafemale['premium2018'])[0, 1]

# 计算 gamma
datafemale['gamma232'] = 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
gamma232= datafemale['gamma232'].mean()
print(gamma232)

970.7524008334406

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
}

# ===== 各异质性的 sigma =====
sigma_map = {
    2: 2.6,  3: 3.4,  4: 3.5,  5: 2.7,
    6: 3.5,  7: 2.7,  8: 3.0,  9: 3.4,
    10: 3.8, 11: 2.8, 12: 3.0, 13: 3.5,
    14: 2.7, 15: 3.0, 16: 3.6, 17: 2.6,
    18: 3.0, 19: 3.6, 20: 2.7, 21: 3.0, 22: 3.5
}

def _to_num(s):  # 小工具：转数值
    return pd.to_numeric(s, errors="coerce")

def compute_gamma_variant2(sub: pd.DataFrame, sigma_val: float) -> float:
    """
    你的公式（可变 σ）：
      E_m0,E_m1 用子样本均值；
      Ec0 = E[c0^(1-σ)]、Ec1 = E[c1^(1-σ)];
      cov0 = Cov( c0^(1-σ)/Ec0, (r0-r1)*m0 + p15 - p18 ), cov1 同理；
      gamma_i = |p15-p18| + |0.5*(r0-r1)*(E_m0+E_m1)| + 0.5*cov0 + 0.5*cov1；
      返回子样本内 gamma_i 的均值。
    """
    if sub.empty:
        return float("nan")

    # 数值化
    for col in ["c0","c1","r0","r1","m0","m1","premium2015","premium2018"]:
        sub[col] = _to_num(sub[col])

    # 子样本均值
    E_m0 = sub["m0"].mean()
    E_m1 = sub["m1"].mean()

    power =  - float(sigma_val)
    c0_pow = sub["c0"] ** power
    c1_pow = sub["c1"] ** power
    Ec0 = c0_pow.mean()
    Ec1 = c1_pow.mean()

    # 协方差
    y0 = (sub["r0"] - sub["r1"]) * sub["m0"] + (sub["premium2015"] - sub["premium2018"])
    y1 = (sub["r0"] - sub["r1"]) * sub["m1"] + (sub["premium2015"] - sub["premium2018"])
    cov0 = np.cov(c0_pow / Ec0, y0)[0, 1]
    cov1 = np.cov(c1_pow / Ec1, y1)[0, 1]

    gamma_series = (sub["premium2015"] - sub["premium2018"]).abs() \
                 + (0.5 * (sub["r0"] - sub["r1"]) * (E_m0 + E_m1)).abs() \
                 + 0.5 * cov0 + 0.5 * cov1

    return float(gamma_series.mean())

# ===== 批量计算：gamma222 … gamma2222 =====
results = {}
for idx in range(2, 23):
    cond_fn = conds[idx]
    # 取子样本（若该条件列不存在，会抛错 -> 该组置为空）
    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 = pd.DataFrame(columns=data.columns)

    name = f"gamma2{idx}2"
    results[name] = compute_gamma_variant2(sub, sigma_map[idx])
    globals()[name] = results[name]  # 可选：注册为同名变量

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

gamma222 = 806.2642603893461
gamma232 = 970.7524008334406
gamma242 = 965.7166817845402
gamma252 = 541.4219608892041
gamma262 = 895.4592545354282
gamma272 = 443.7762743250522
gamma282 = 638.479504002872
gamma292 = 933.2970096494091
gamma2102 = 936.2440886871743
gamma2112 = 1035.3304448109861
gamma2122 = 994.600066804255
gamma2132 = 649.6175048828753
gamma2142 = 850.429809723793
gamma2152 = 974.6068796763933
gamma2162 = 793.0948494149856
gamma2172 = 849.3248136429944
gamma2182 = 584.6022986316312
gamma2192 = 1004.9526781209704
gamma2202 = 1096.1588816011183
gamma2212 = 746.2723252894853
gamma2222 = 857.3239372490277

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
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
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

#计算dh/dm 15
import pandas as pd
import statsmodels.api as sm
e1 = pd.merge(c0m0,hsf15,on="ID",how="inner")
e1= e1[['m0','hsf15']].copy()
# 删除包含 NaN 或 inf 的行
e1= e1.replace([np.inf, -np.inf], np.nan).dropna()
# 删除包含 0 的行
e1 = e1[(e1['hsf15']!=0) & (e1['m0']!=0)]
# 自变量（X）和因变量（Y）
X = e1['m0']
Y = e1['hsf15']

# 在 X 中添加常数项，以便进行 OLS 回归
X = sm.add_constant(X)

# 拟合 OLS 回归模型
model = sm.OLS(Y, X).fit()

# 输出回归结果
print(model.summary())

# 提取回归系数
coefficients = model.params

# 保存特定自变量的回归系数
h_m15 = coefficients['m0'] 
h_m15

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  hsf15   R-squared:                       0.000
Model:                            OLS   Adj. R-squared:                  0.000
Method:                 Least Squares   F-statistic:                     1.297
Date:                Mon, 29 Dec 2025   Prob (F-statistic):              0.255
Time:                        21:23:42   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.

np.float64(-1.7062284276591698e-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.579
Date:                Mon, 29 Dec 2025   Prob (F-statistic):              0.209
Time:                        21:23:42   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.

np.float64(-1.3590936773055398e-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^(-sigma)) 和 E(c1^(-sigma))
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['gamma213'] = 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
gamma213 = dataa['gamma213'].mean()
print(gamma213)
dataa

689.4542339864197

	ID	c0	m0	householdID	communityID	c1	m1	gender	age	marriage	...	r0	r1	r0adjust	r1adjust	policyintergration2015	policyintergration2018	district	GDPgrowthrate	urban_nbs	gamma213
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

#健康based异质性—男性
dataamale = dataa.loc[dataa['gender'] == 1].copy()
sigamamale = 2.6
import pandas as pd
import numpy as np

# 计算 E(m0) 和 E(m1)
E_m0 = dataamale['m0'].mean()
E_m1 = dataamale['m1'].mean()

# 计算 E(c0^(-sigma)) 和 E(c1^(-sigma))
E_c0_inv2 = (dataamale['c0']**(-sigamamale)).mean()
E_c1_inv2 = (dataamale['c1']**(-sigamamale)).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['gamma223'] = 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
gamma223= dataamale['gamma223'].mean()
print(gamma223)

679.9212294518087

#健康based异质性—女性
dataafemale = dataa.loc[dataa['gender'] == 0].copy()
sigamafemale = 3.4
import pandas as pd
import numpy as np

# 计算 E(m0) 和 E(m1)
E_m0 = dataafemale['m0'].mean()
E_m1 = dataafemale['m1'].mean()

# 计算 E(c0^(-sigma)) 和 E(c1^(-sigma))
E_c0_inv2 = (dataafemale['c0']**(-sigamafemale)).mean()
E_c1_inv2 = (dataafemale['c1']**(-sigamafemale)).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['gamma233'] = 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
gamma233 = dataafemale['gamma233'].mean()
print(gamma233)

772.8992376511646

import numpy as np
import pandas as pd

PHI = 0.019743  # phi_tilde

# ===== 异质性条件 =====
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
}

# ===== 每组的 σ =====
sigma_map = {
    2: 2.6,  3: 3.4,  4: 3.5,  5: 2.7,
    6: 3.5,  7: 2.7,  8: 3.0,  9: 3.4,
    10: 3.8, 11: 2.8, 12: 3.0, 13: 3.5,
    14: 2.7, 15: 3.0, 16: 3.6, 17: 2.6,
    18: 3.0, 19: 3.6, 20: 2.7, 21: 3.0, 22: 3.5
}

# ===== 小工具 =====
def _to_num(s):  # 数值化
    return pd.to_numeric(s, errors="coerce")

def _safe_cov(x, y):  # 样本协方差，自动清理 NaN/Inf
    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 与子样本对齐：
    - 若 h_full 是 Series 且 index == dataa_full.index：按 index 对齐
    - 若 h_full 的 index 是 ID：按 ID reindex
    - 若是数组：包成 Series(index=dataa_full.index) 后再切
    """
    if isinstance(h_full, pd.Series):
        if h_full.index.equals(dataa_full.index):
            return h_full.loc[data_sub.index]
        if id_col in data_sub.columns and h_full.index.isin(data_sub[id_col]).any():
            s = h_full.reindex(data_sub[id_col])
            s.index = data_sub.index
            return s
        try:
            return h_full.loc[data_sub.index]
        except Exception:
            return pd.Series(np.nan, index=data_sub.index)
    # array-like
    try:
        base = pd.Series(h_full, index=dataa_full.index)
        return base.loc[data_sub.index]
    except Exception:
        return pd.Series(np.nan, index=data_sub.index)

def compute_gamma_health_sigma(dataa_full: pd.DataFrame, data_sub: pd.DataFrame, sigma_val: float,
                               h_m15, h_m18) -> float:
    """
    对某一子样本按给定 σ 计算：
      Ec0 = E[c0^(1-σ)]，Ec1 = E[c1^(1-σ)]
      cov0 = Cov( (PHI*h15)/(Ec0*r0), (r0-r1)*m0 + p15 - p18 )
      cov1 = Cov( (PHI*h18)/(Ec1*r1), (r0-r1)*m1 + p15 - p18 )
      γ = mean( |p15-p18| + |0.5*(r0-r1)*(E[m0]+E[m1])| ) + 0.5*(cov0+cov1)
    """
    if data_sub.empty:
        return float("nan")

    # 数值化
    for col in ["c0","c1","r0","r1","m0","m1","premium2015","premium2018"]:
        data_sub[col] = _to_num(data_sub[col])

    # 计算 E[c^(-σ)]
    power = - float(sigma_val)
    # 为避免 0 的负幂 -> inf，先把 0 当作缺失
    c0_pow = (data_sub["c0"].replace(0, np.nan)) ** power
    c1_pow = (data_sub["c1"].replace(0, np.nan)) ** power
    Ec0 = c0_pow.replace([np.inf, -np.inf], np.nan).mean()
    Ec1 = c1_pow.replace([np.inf, -np.inf], np.nan).mean()

    # 对齐健康项
    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"]

    # 两个协方差
    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)

    # E[m0], E[m1]
    Em0 = m0.mean()
    Em1 = m1.mean()

    gamma_series = (p15 - p18).abs() + (0.5 * (data_sub["r0"] - data_sub["r1"]) * (Em0 + Em1)).abs()
    gamma_val = float(gamma_series.mean() + 0.5*cov0 + 0.5*cov1)
    return gamma_val

# ===== 批量计算：gamma223 … gamma2223 =====
_results = {}
for idx in range(2, 23):
    cond_fn = conds[idx]
    # 取子样本；若条件列缺失则该组返回 NaN（不影响其它组）
    try:
        mask = cond_fn(dataa)
        sub = dataa.loc[mask].copy() if isinstance(mask, pd.Series) and len(mask)==len(dataa) else pd.DataFrame(columns=dataa.columns)
    except Exception:
        sub = pd.DataFrame(columns=dataa.columns)

    name = f"gamma2{idx}3"
    _results[name] = compute_gamma_health_sigma(dataa_full=dataa, data_sub=sub,
                                                sigma_val=sigma_map[idx],
                                                h_m15=h_m15, h_m18=h_m18)
    globals()[name] = _results[name]  # 可选：注册为同名变量

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

gamma223 = 679.9212294518087
gamma233 = 772.8992376511646
gamma243 = 964.4520934881297
gamma253 = 482.75095269985326
gamma263 = 870.5672626997266
gamma273 = 362.56619289490567
gamma283 = 671.0503320289147
gamma293 = 755.4755004388335
gamma2103 = 1392.6896145700953
gamma2113 = 597.4245950966861
gamma2123 = 763.740196317552
gamma2133 = 729.5221044951721
gamma2143 = 617.6144223711593
gamma2153 = 795.3126953827027
gamma2163 = 873.2308196485699
gamma2173 = 605.2717675894386
gamma2183 = 555.6822076078611
gamma2193 = 1009.6456471211893
gamma2203 = 681.1315362913325
gamma2213 = 813.6137702363887
gamma2223 = 791.10194346691

用完全信息法求解

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
...	...	...	...	...	...	...	...
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

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

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

1996.3867642406108

	ID	c0	m0	c1	m1	hsf15	hsf18	B_bar	gamma211
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

#异质性男性
import numpy as np
import pandas as pd
d4 = pd.merge(d1, hsf18full, on="ID", how="inner")
dmale = d4.loc[d4['gender'] == 1].copy()

# 参数
sigmamale = 2.6
phi_tilde = 0.019743

dmale["B_bar"] = (dmale["c0"]**(1 - sigmamale)) + (1 - sigmamale) * phi_tilde * (dmale["hsf15"] - dmale["hsf18"])
dmale["gamma221"] = dmale["c1"] - dmale["B_bar"]**(1 / (1 - sigmamale))

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

1968.0150747788323

#异质性女性
import numpy as np
import pandas as pd
d4 = pd.merge(d1, hsf18full, on="ID", how="inner")
dfemale = d4.loc[d4['gender'] == 0].copy()

# 参数
sigmafemale = 3.4
phi_tilde = 0.019743

dfemale["B_bar"] = (dfemale["c0"]**(1 - sigmafemale)) + (1 - sigmafemale) * phi_tilde * (dfemale["hsf15"] - dfemale["hsf18"])
dfemale["gamma231"] = dfemale["c1"] - dfemale["B_bar"]**(1 / (1 - sigmafemale))

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

2030.7369827386915

import numpy as np
import pandas as pd

# 固定合表
d4 = pd.merge(d1, hsf18full, on="ID", how="inner")

# 参与计算的列转数值（若本来就是数值类型，这步不会改变结果）
for col in ["c0", "c1", "hsf15", "hsf18"]:
    if col in d4.columns:
        d4[col] = pd.to_numeric(d4[col], errors="coerce")

# —— 异质性条件 ——
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
}

# —— 每组的 sigma ——
sigma_map = {
    2: 2.6,  3: 3.4,  4: 3.5,  5: 2.7,
    6: 3.5,  7: 2.7,  8: 3.0,  9: 3.4,
    10: 3.8, 11: 2.8, 12: 3.0, 13: 3.5,
    14: 2.7, 15: 3.0, 16: 3.6, 17: 2.6,
    18: 3.0, 19: 3.6, 20: 2.7, 21: 3.0, 22: 3.5
}

phi_tilde = 0.019743

def compute_gamma_group(df_sub: pd.DataFrame, sigma_val: float) -> float:
    """
    子样本逐行计算：
      B_i = c0^(1-σ) + (1-σ)*phi*(hsf15 - hsf18)
      gamma_i = c1 - B_i^(1/(1-σ))
    返回子样本内 gamma_i 的均值（忽略 NaN）
    """
    if df_sub.empty:
        return float("nan")
    power = 1.0 - float(sigma_val)
    inv_power = 1.0 / power

    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(skipna=True))

# —— 批量：gamma221 … gamma2221 ——
results = {}
for idx in range(2, 23):
    # 取子样本；若条件列缺失则该组置空
    try:
        mask = conds[idx](d4)
        sub = d4.loc[mask].copy() if isinstance(mask, pd.Series) and len(mask)==len(d4) else pd.DataFrame(columns=d4.columns)
    except Exception:
        sub = pd.DataFrame(columns=d4.columns)

    name = f"gamma2{idx}1"
    results[name] = compute_gamma_group(sub, sigma_map[idx])
    globals()[name] = results[name]  # 可选：注册为同名变量

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

gamma221 = 1968.0150747788323
gamma231 = 2030.7369827386915
gamma241 = 2091.9530666092396
gamma251 = 1516.3798526359226
gamma261 = 2009.9144694593858
gamma271 = 751.9705731491155
gamma281 = 2513.485763337119
gamma291 = 1492.4551887750622
gamma2101 = 964.8652631218005
gamma2111 = 2318.2529510457402
gamma2121 = 2056.9974605815864
gamma2131 = 1614.9215720119535
gamma2141 = 2312.557965506386
gamma2151 = 1938.9117698957555
gamma2161 = 1560.412555197337
gamma2171 = 2566.3049550302685
gamma2181 = 2138.1528481056553
gamma2191 = 1882.7784764569944
gamma2201 = 2677.534672903796
gamma2211 = 2276.608300019149
gamma2221 = 1740.8856369742805

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

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

data = [
[gamma211, gamma212, gamma213],
[gamma221, gamma222, gamma223],
[gamma231, gamma232, gamma233],
[gamma241, gamma242, gamma243],
[gamma251, gamma252, gamma253],
[gamma261, gamma262, gamma263],
[gamma271, gamma272, gamma273],
[gamma281, gamma282, gamma283],
[gamma291, gamma292, gamma293],
[gamma2101, gamma2102, gamma2103],
[gamma2111, gamma2112, gamma2113],
[gamma2121, gamma2122, gamma2123],
[gamma2131, gamma2132, gamma2133],
[gamma2141, gamma2142, gamma2143],
[gamma2151, gamma2152, gamma2153],
[gamma2161, gamma2162, gamma2163],
[gamma2171, gamma2172, gamma2173],
[gamma2181, gamma2182, gamma2183],
[gamma2191, gamma2192, gamma2193],
[gamma2201, gamma2202, gamma2203],
[gamma2211, gamma2212, gamma2213],
[gamma2221, gamma2222, gamma2223],
]

# 2) 多级列索引
cols = pd.MultiIndex.from_tuples([
    ("完全信息方法",""),
    ("最优化方法","仅假设效用函数\n的消费部分"),
    ("最优化方法","仅假设效用函数\n的健康部分"),
])

df = pd.DataFrame(data, index=rows, columns=cols)

# 3) 分组起始行（加粗横线）
group_starts = {
    "男性",           # 性别组
    "有配偶",         # 婚姻组
    "有子女",         # 子女组
    "45 岁—59 岁",    # 年龄组
    "东部",           # 地区组
    "健康状况较好",   # 健康组
    "较高收入",       # 收入组
    "教育程度较高"    # 教育组
}

def row_borders(row):
    label = row.name
    if label in group_starts:
        return ['border-top: 2px solid #4a4a4a'] * len(row)
    return [''] * len(row)

# 4) 样式与展示
styler = (
    df.style
      .set_table_styles([
          {'selector': 'th.col_heading.level0',
           'props': [('font-weight', '700'),
                     ('border-bottom','1px solid #4a4a4a')]},
          {'selector': 'th.col_heading.level1',
           'props': [('font-weight', '700')]},
          {'selector': 'th.row_heading',
           'props': [('font-weight', '700')]},
          {'selector': 'table',
           'props': [('border-collapse','collapse'),
                     ('font-family','-apple-system,BlinkMacSystemFont,Segoe UI,Roboto,PingFang SC,Helvetica,Arial')]}
      ])
      .format(precision=0)
      .set_properties(**{
          'text-align': 'center',
          'padding': '6px',
          'border':'1px solid #a0a0a0'
      })
      .apply(row_borders, axis=1)
)

# 在 Jupyter 中显示
styler

# （可选）导出为 HTML 文件
# with open("表格_完全信息与最优化方法.html", "w", encoding="utf-8") as f:
#     f.write(styler.to_html())

	完全信息方法	最优化方法
		仅假设效用函数 的消费部分	仅假设效用函数 的健康部分
全样本	1996	894	689
男性	1968	806	680
女性	2031	971	773
有配偶	2092	966	964
无配偶	1516	541	483
有子女	2010	895	871
无子女	752	444	363
小于59 岁	2513	638	671
60 岁—79 岁	1492	933	755
80 岁及以上	965	936	1393
东部	2318	1035	597
中部	2057	995	764
西部	1615	650	730
健康状况较好	2313	850	618
健康状况中等	1939	975	795
健康状况较差	1560	793	873
较高收入	2566	849	605
中等收入	2138	585	556
较低收入	1883	1005	1010
教育程度较高	2678	1096	681
教育程度中等	2277	746	814
教育程度较低	1741	857	791