导入数据

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(1/c0) 和 E(1/c1)
E_c0_inv2 = (1/data['c0']).mean()
E_c1_inv2 = (1/data['c1']).mean()

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

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

841.987535959292

	ID	c0	m0	householdID	communityID	c1	m1	gender	age	marriage	...	r1	r0adjust	r1adjust	policyintergration2015	policyintergration2018	district	GDPgrowthrate	urban_nbs	gamma612	gamma312
0	64033321002	112.216	60.0	640333210	640333	123.670	0.0	0.0	59.0	1.0	...	0.700	0.544192	0.660903	0.0	1.0	east	0.118005	Rural	652.775511	652.775511
1	64033327002	1029.200	6000.0	640333270	640333	1054.764	21000.0	0.0	62.0	1.0	...	0.700	0.544192	0.660903	0.0	1.0	east	0.118005	Rural	652.775511	652.775511
2	64033325001	1062.400	3000.0	640333250	640333	78.020	100.0	1.0	66.0	1.0	...	0.700	0.544192	0.660903	0.0	1.0	east	0.118005	Rural	652.775511	652.775511
3	64033322001	592.620	2000.0	640333220	640333	859.050	17050.0	1.0	63.0	1.0	...	0.700	0.544192	0.660903	0.0	1.0	east	0.118005	Rural	652.775511	652.775511
4	64033330002	2058.400	2000.0	640333300	640333	4025.500	2000.0	1.0	59.0	1.0	...	0.700	0.544192	0.660903	0.0	1.0	east	0.118005	Rural	652.775511	652.775511
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
3882	89676104001	2315.700	840.0	896761040	896761	891.420	8000.0	1.0	61.0	1.0	...	0.725	0.624462	0.685108	0.0	1.0	west	0.284050	Rural	701.272074	701.272074
3883	89676114002	1935.145	300.0	896761140	896761	49.800	0.0	0.0	56.0	0.0	...	0.725	0.624462	0.685108	0.0	1.0	west	0.284050	Rural	701.272074	701.272074
3884	89676118001	2466.096	1000.0	896761180	896761	661.095	3000.0	0.0	73.0	0.0	...	0.725	0.624462	0.685108	0.0	1.0	west	0.284050	Rural	701.272074	701.272074
3885	89676115001	10721.940	800.0	896761150	896761	11638.260	2000.0	1.0	55.0	1.0	...	0.725	0.624462	0.685108	0.0	1.0	west	0.284050	Rural	701.272074	701.272074
3886	89676124001	268.422	500.0	896761240	896761	313.242	101.0	1.0	69.0	0.0	...	0.725	0.624462	0.685108	0.0	1.0	west	0.284050	Rural	701.272074	701.272074

3887 rows × 28 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^(-2)) 和 E(c1^(-2))
E_c0_inv2 = (1/datamale['c0']).mean()
E_c1_inv2 = (1/datamale['c1']).mean()

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

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

802.7906600812382

#异质性女性 混合截面
#最优化方法——消费计算 
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^(-2)) 和 E(c1^(-2))
E_c0_inv2 = (1/datafemale['c0']).mean()
E_c1_inv2 = (1/datafemale['c1']).mean()

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

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

881.1201345894027

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

def compute_gamma_inv1(sub: pd.DataFrame) -> float:
    """
    本方法（对应 gamma3?2）：
      Ec0 = E[1/c0]，Ec1 = E[1/c1]
      cov0 = Cov( (1/c0)/Ec0 , (r0-r1)*m0 + p15 - p18 )
      cov1 = Cov( (1/c1)/Ec1 , (r0-r1)*m1 + p15 - p18 )
      γ_i = |p15-p18| + |0.5*(r0-r1)*(E[m0]+E[m1])| + 0.5*cov0 + 0.5*cov1
      返回子样本内 γ_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()

    # 1/c，避免除以 0：把 0 当缺失
    c0_inv = 1.0 / sub["c0"].replace(0, np.nan)
    c1_inv = 1.0 / sub["c1"].replace(0, np.nan)
    Ec0 = c0_inv.replace([np.inf, -np.inf], np.nan).mean()
    Ec1 = c1_inv.replace([np.inf, -np.inf], np.nan).mean()

    y0 = (sub["r0"] - sub["r1"]) * sub["m0"] + (sub["premium2015"] - sub["premium2018"])
    y1 = (sub["r0"] - sub["r1"]) * sub["m1"] + (sub["premium2015"] - sub["premium2018"])

    cov0 = _safe_cov(c0_inv / Ec0, y0) if pd.notna(Ec0) and Ec0 != 0 else 0.0
    cov1 = _safe_cov(c1_inv / Ec1, y1) if pd.notna(Ec1) and Ec1 != 0 else 0.0

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

# —— 批量计算：gamma322 … gamma3222 ——
results = {}
for idx in range(2, 23):
    try:
        mask = conds[idx](data)
        sub = data.loc[mask].copy() if isinstance(mask, pd.Series) and len(mask)==len(data) else pd.DataFrame(columns=data.columns)
    except Exception:
        sub = pd.DataFrame(columns=data.columns)

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

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

gamma322 = 802.7906600812382
gamma332 = 881.1201345894027
gamma342 = 887.7108284615841
gamma352 = 576.663420801471
gamma362 = 846.9163570623545
gamma372 = 394.1618493492495
gamma382 = 740.2572102706334
gamma392 = 878.2989375956844
gamma3102 = 771.6318089577574
gamma3112 = 832.6376734492181
gamma3122 = 954.152803417741
gamma3132 = 658.8869816345907
gamma3142 = 787.623935492869
gamma3152 = 966.6883145675897
gamma3162 = 695.2482563497338
gamma3172 = 749.2273524859135
gamma3182 = 614.4296042994745
gamma3192 = 934.3856349227033
gamma3202 = 810.1788592321084
gamma3212 = 923.318782424085
gamma3222 = 808.3482289017528

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, 13 Oct 2025   Prob (F-statistic):              0.255
Time:                        01:34:37   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, 13 Oct 2025   Prob (F-statistic):              0.209
Time:                        01:34:37   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(1/c0) 和 E(1/c1)
E_c0_inv2 = (1/dataa['c0']).mean()
E_c1_inv2 = (1/dataa['c1']).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['gamma313'] = 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
gamma313 = dataa['gamma313'].mean()
print(gamma313)
dataa

665.7478215297924

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

3887 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_inv2 = (1/dataamale['c0']).mean()
E_c1_inv2 = (1/dataamale['c1']).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['gamma323'] = 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
gamma323 = dataamale['gamma323'].mean()
print(gamma323)

675.5628844320822

#健康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_inv2 = (1/dataafemale['c0']).mean()
E_c1_inv2 = (1/dataafemale['c1']).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['gamma333'] = 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
gamma333 = dataafemale['gamma333'].mean()
print(gamma333)

653.9915283097725

import numpy as np
import pandas as pd

PHI = 0.019743  # 0.019743 常数

# —— 异质性条件——
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):
    """丢弃 NaN/Inf 后算样本协方差（与你 np.cov 默认口径一致）"""
    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])[0,1])
    return 0.0

def _align_health_for_subset(dataa_full: pd.DataFrame, data_sub: pd.DataFrame, h_full, id_col="ID"):
    """
    将 h_m15 / h_m18 与子样本按 index 或 ID 对齐；h_full 可为 Series 或 array-like。
    """
    if isinstance(h_full, pd.Series):
        # 与 dataa_full 同 index
        if h_full.index.equals(dataa_full.index):
            return h_full.loc[data_sub.index]
        # 若 h_full 的索引可能是 ID
        if id_col in data_sub.columns and h_full.index.isin(data_sub[id_col]).any():
            s = h_full.reindex(data_sub[id_col])
            s.index = data_sub.index
            return s
        # 兜底按 index 选
        try:
            return h_full.loc[data_sub.index]
        except Exception:
            return pd.Series(np.nan, index=data_sub.index)
    # array-like：包成与 dataa_full 对齐的 Series 再切
    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_inv1(dataa_full: pd.DataFrame, data_sub: pd.DataFrame, h_m15, h_m18) -> float:
    """
    本方法（对应 gamma3?3）：
      Ec0 = E[1/c0]，Ec1 = E[1/c1]
      cov0 = Cov( (PHI*h15)/(Ec0*r0) , (r0-r1)*m0 + p15 - p18 )
      cov1 = Cov( (PHI*h18)/(Ec1*r1) , (r0-r1)*m1 + p15 - p18 )
      γ_i  = |p15-p18| + |0.5*(r0-r1)*(E[m0]+E[m1])| + 0.5*cov0 + 0.5*cov1
      返回子样本内 γ_i 的均值
    """
    if data_sub.empty:
        return float("nan")

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

    # E[m0], E[m1]
    E_m0 = data_sub["m0"].mean()
    E_m1 = data_sub["m1"].mean()

    # E[1/c]（把无穷替换为 NaN 再求均值，避免 1/0 影响）
    c0_inv = 1.0 / data_sub["c0"]
    c1_inv = 1.0 / data_sub["c1"]
    Ec0 = c0_inv.replace([np.inf, -np.inf], np.nan).mean()
    Ec1 = c1_inv.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)

    # y0,y1
    y0 = (data_sub["r0"] - data_sub["r1"]) * data_sub["m0"] + (data_sub["premium2015"] - data_sub["premium2018"])
    y1 = (data_sub["r0"] - data_sub["r1"]) * data_sub["m1"] + (data_sub["premium2015"] - data_sub["premium2018"])

    # 两个协方差
    r0 = data_sub["r0"]
    r1 = data_sub["r1"]
    cov0 = _safe_cov((PHI * h15) / (Ec0 * r0), y0) if pd.notna(Ec0) and Ec0 != 0 else 0.0
    cov1 = _safe_cov((PHI * h18) / (Ec1 * r1), y1) if pd.notna(Ec1) and Ec1 != 0 else 0.0

    gamma_series = (data_sub["premium2015"] - data_sub["premium2018"]).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())

# —— 批量计算：gamma323 … gamma3223 ——
results = {}
for idx in range(2, 23):
    try:
        mask = conds[idx](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"gamma3{idx}3"
    results[name] = compute_gamma_health_inv1(dataa_full=dataa, data_sub=sub, h_m15=h_m15, h_m18=h_m18)
    globals()[name] = results[name]  # 可选：注册为同名变量

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

gamma323 = 675.5628844320822
gamma333 = 653.9915283097725
gamma343 = 711.7527340731003
gamma353 = 478.53226004335176
gamma363 = 669.9089945050528
gamma373 = 364.73161872871424
gamma383 = 646.919021748351
gamma393 = 679.7402110445717
gamma3103 = 539.7100207150593
gamma3113 = 594.8823272349582
gamma3123 = 741.652210399491
gamma3133 = 584.1947275016881
gamma3143 = 616.7061952195254
gamma3153 = 763.3067968734329
gamma3163 = 560.3194837995072
gamma3173 = 589.0182650573405
gamma3183 = 540.942205343509
gamma3193 = 715.2644872566253
gamma3203 = 650.674324132809
gamma3213 = 748.8317568381137
gamma3223 = 637.7790899630422

用完全信息法求解

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
# 参数
PHI = 0.019743
d3["gamma311"] = d3["c1"] * (PHI * (d3["hsf18"] - d3["hsf15"]) - np.log(d3["c0"] / d3["c1"]))

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

2069.455364572046

	ID	c0	m0	c1	m1	hsf15	hsf18	gamma311
0	64033321002	112.216	60.0	123.670	0.0	54.436016	10.551840	-95.128712
1	64033327002	1029.200	6000.0	1054.764	21000.0	28.096293	39.754467	268.651118
2	64033322001	592.620	2000.0	859.050	17050.0	21.303569	59.803845	971.916021
3	64033330002	2058.400	2000.0	4025.500	2000.0	34.523055	68.626626	5410.380456
4	64033341001	5436.500	500.0	1806.578	2500.0	50.384267	78.657775	-981.870016
...	...	...	...	...	...	...	...	...
3750	89676104001	2315.700	840.0	891.420	8000.0	27.522919	42.058152	-595.185530
3751	89676114002	1935.145	300.0	49.800	0.0	35.235805	15.917436	-201.257990
3752	89676118001	2466.096	1000.0	661.095	3000.0	32.251424	54.653869	-577.930962
3753	89676115001	10721.940	800.0	11638.260	2000.0	32.757347	73.821948	10389.989777
3754	89676124001	268.422	500.0	313.242	101.0	34.706036	95.825352	426.351916

3755 rows × 8 columns

#异质性计算男性
import numpy as np
import pandas as pd
# 参数
PHI = 0.019743

d4 = pd.merge(d1, hsf18full, on="ID", how="inner")
dmale = d4.loc[d4['gender'] == 1].copy()

dmale["gamma321"] = dmale["c1"] * (PHI * (dmale["hsf18"] - dmale["hsf15"]) - np.log(dmale["c0"] / dmale["c1"]))

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

1972.9361597607306

#异质性计算女性
import numpy as np
import pandas as pd
# 参数
PHI = 0.019743

d4 = pd.merge(d1, hsf18full, on="ID", how="inner")
dfemale = d4.loc[d4['gender'] == 0].copy()

dfemale["gamma331"] = dfemale["c1"] * (PHI * (dfemale["hsf18"] - dfemale["hsf15"]) - np.log(dfemale["c0"] / dfemale["c1"]))

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

2172.511938431705

import numpy as np
import pandas as pd

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

PHI = 0.019743

# —— 异质性条件——
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 compute_gamma(sub: pd.DataFrame) -> float:
    """
    γ_i = c1 * ( PHI*(hsf18 - hsf15) - ln(c0/c1) )
    返回子样本内 γ_i 的均值
    """
    if sub.empty:
        return float("nan")
    # 若需要更稳健可数值化；为尽量贴近你原始结果，这里不强制转型
    gamma_series = sub["c1"] * (PHI * (sub["hsf18"] - sub["hsf15"]) - np.log(sub["c0"] / sub["c1"]))
    return float(gamma_series.mean())

# —— 批量计算：gamma321 … gamma3221 ——
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"gamma3{idx}1"
    results[name] = compute_gamma(sub)
    globals()[name] = results[name]  # 可选：注册为同名变量，便于后续直接用

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

gamma321 = 1972.9361597607306
gamma331 = 2172.511938431705
gamma341 = 2133.875674113165
gamma351 = 1810.9140421680336
gamma361 = 2089.855534379301
gamma371 = 579.8554327323171
gamma381 = 2722.573800517208
gamma391 = 1519.6215845307818
gamma3101 = 816.5738606477122
gamma3111 = 2828.1348185073775
gamma3121 = 2138.774183188456
gamma3131 = 1400.2565528289158
gamma3141 = 2362.867406843136
gamma3151 = 2091.3356304820554
gamma3161 = 1631.1663360643925
gamma3171 = 2620.5894489874076
gamma3181 = 2070.3434949871958
gamma3191 = 2016.9799313921617
gamma3201 = 3566.7253958692945
gamma3211 = 2454.8369532172455
gamma3221 = 1720.6506749898997

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

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

data = [
[gamma311, gamma312, gamma313],
[gamma321, gamma322, gamma323],
[gamma331, gamma332, gamma333],
[gamma341, gamma342, gamma343],
[gamma351, gamma352, gamma353],
[gamma361, gamma362, gamma363],
[gamma371, gamma372, gamma373],
[gamma381, gamma382, gamma383],
[gamma391, gamma392, gamma393],
[gamma3101, gamma3102, gamma3103],
[gamma3111, gamma3112, gamma3113],
[gamma3121, gamma3122, gamma3123],
[gamma3131, gamma3132, gamma3133],
[gamma3141, gamma3142, gamma3143],
[gamma3151, gamma3152, gamma3153],
[gamma3161, gamma3162, gamma3163],
[gamma3171, gamma3172, gamma3173],
[gamma3181, gamma3182, gamma3183],
[gamma3191, gamma3192, gamma3193],
[gamma3201, gamma3202, gamma3203],
[gamma3211, gamma3212, gamma3213],
[gamma3221, gamma3222, gamma3223],
]

# 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

	完全信息方法	最优化方法
		仅假设效用函数 的消费部分	仅假设效用函数 的健康部分
全样本	2069	842	666
男性	1973	803	676
女性	2173	881	654
有配偶	2134	888	712
无配偶	1811	577	479
有子女	2090	847	670
无子女	580	394	365
小于59 岁	2723	740	647
60 岁—79 岁	1520	878	680
80 岁及以上	817	772	540
东部	2828	833	595
中部	2139	954	742
西部	1400	659	584
健康状况较好	2363	788	617
健康状况中等	2091	967	763
健康状况较差	1631	695	560
较高收入	2621	749	589
中等收入	2070	614	541
较低收入	2017	934	715
教育程度较高	3567	810	651
教育程度中等	2455	923	749
教育程度较低	1721	808	638