导入数据¶

In [2]:
import pandas as pd
c0m0full = pd.read_csv(r'D:\论文\最后一波一鼓作气\数据\c0m0full.csv')  
c1m1full = pd.read_csv(r'D:\论文\最后一波一鼓作气\数据\c1m1full.csv')  
hsf15full = pd.read_csv(r'D:\论文\最后一波一鼓作气\数据\hsf15full.csv')  
hsf18full = pd.read_csv(r'D:\论文\最后一波一鼓作气\数据\hsf18full.csv')  
city_gender_age_premium_ratio_district15 = pd.read_csv(r'D:\论文\最后一波一鼓作气\数据\city_gender_age_premium_ratio_district15.csv')
In [3]:
#删除城市样本
c1m1full = c1m1full[c1m1full['urban_nbs'] != 'Urban']
#确保15年未整合，18年整合了
c1m1full = c1m1full[(c1m1full['policyintergration2015']==0.0) & (c1m1full['policyintergration2018']==1.0)]
c1m1= c1m1full[['ID', 'c1','m1']]

#删除城市样本
c0m0full = c0m0full[c0m0full['urban_nbs'] != 'Urban']
#确保15年未整合，18年整合了
c0m0full = c0m0full[(c0m0full['policyintergration2015']==0.0) & (c0m0full['policyintergration2018']==1.0)]
c0m0= c0m0full[['ID', 'c0','m0']]

#删除城市样本
hsf15full = hsf15full[hsf15full['urban_nbs'] != 'Urban']
#确保15年未整合，18年整合了
hsf15full = hsf15full[(hsf15full['policyintergration2015']==0.0) & (hsf15full['policyintergration2018']==1.0)]
hsf15= hsf15full[['ID', 'hsf15']]

#删除城市样本
hsf18full = hsf18full[hsf18full['urban_nbs'] != 'Urban']
#确保15年未整合，18年整合了
hsf18full = hsf18full[(hsf18full['policyintergration2015']==0.0) & (hsf18full['policyintergration2018']==1.0)]
hsf18= hsf18full[['ID', 'hsf18']]

最优化方法¶

consumption-based：合并数据¶

In [4]:
data = pd.merge(c0m0, c1m1full,on="ID", how="inner")
data
Out[4]:
ID c0 m0 householdID communityID c1 m1 gender age marriage ... premium2018 r0 r1 r0adjust r1adjust policyintergration2015 policyintergration2018 district GDPgrowthrate urban_nbs
0 64033321002 112.216 60.0 640333210 640333 123.670 0.0 0.0 59.0 1.0 ... 180.0 0.6167 0.700 0.544192 0.660903 0.0 1.0 east 0.118005 Rural
1 64033327002 1029.200 6000.0 640333270 640333 1054.764 21000.0 0.0 62.0 1.0 ... 180.0 0.6167 0.700 0.544192 0.660903 0.0 1.0 east 0.118005 Rural
2 64033325001 1062.400 3000.0 640333250 640333 78.020 100.0 1.0 66.0 1.0 ... 180.0 0.6167 0.700 0.544192 0.660903 0.0 1.0 east 0.118005 Rural
3 64033322001 592.620 2000.0 640333220 640333 859.050 17050.0 1.0 63.0 1.0 ... 180.0 0.6167 0.700 0.544192 0.660903 0.0 1.0 east 0.118005 Rural
4 64033330002 2058.400 2000.0 640333300 640333 4025.500 2000.0 1.0 59.0 1.0 ... 180.0 0.6167 0.700 0.544192 0.660903 0.0 1.0 east 0.118005 Rural
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
3882 89676104001 2315.700 840.0 896761040 896761 891.420 8000.0 1.0 61.0 1.0 ... 220.0 0.6500 0.725 0.624462 0.685108 0.0 1.0 west 0.284050 Rural
3883 89676114002 1935.145 300.0 896761140 896761 49.800 0.0 0.0 56.0 0.0 ... 220.0 0.6500 0.725 0.624462 0.685108 0.0 1.0 west 0.284050 Rural
3884 89676118001 2466.096 1000.0 896761180 896761 661.095 3000.0 0.0 73.0 0.0 ... 220.0 0.6500 0.725 0.624462 0.685108 0.0 1.0 west 0.284050 Rural
3885 89676115001 10721.940 800.0 896761150 896761 11638.260 2000.0 1.0 55.0 1.0 ... 220.0 0.6500 0.725 0.624462 0.685108 0.0 1.0 west 0.284050 Rural
3886 89676124001 268.422 500.0 896761240 896761 313.242 101.0 1.0 69.0 0.0 ... 220.0 0.6500 0.725 0.624462 0.685108 0.0 1.0 west 0.284050 Rural

3887 rows × 26 columns

In [5]:
#最优化方法——消费计算
import pandas as pd
import numpy as np

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

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

#最优化方法——消费计算（表2consumption based混合截面）
gamma312=abs(city_gender_age_premium_ratio_district15['premium2015'].mean() - city_gender_age_premium_ratio_district15['premium2018'].mean()) + abs(0.5 * (city_gender_age_premium_ratio_district15['r0'].mean() - city_gender_age_premium_ratio_district15['r1'].mean()) * (c0m0['m0'].mean() + c1m1['m1'].mean())) + 0.5 * cov_c0_inv2 + 0.5 * cov_c1_inv2
float(gamma312)
Out[5]:
990.1055666990669
In [6]:
#最优化方法——消费计算 平衡面板
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['gamma322'] = 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
gamma322= data['gamma322'].mean()
print(gamma322)
data

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

3887 rows × 27 columns

In [17]:
#最优化方法——消费计算 使用调整之后的报销比例混合截面
import pandas as pd
import numpy as np

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

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

#最优化方法——消费计算（表2consumption based混合截面）
gamma412=abs(city_gender_age_premium_ratio_district15['premium2015'].mean() - city_gender_age_premium_ratio_district15['premium2018'].mean()) + abs(0.5 * (city_gender_age_premium_ratio_district15['r0adjust'].mean() - city_gender_age_premium_ratio_district15['r1adjust'].mean()) * (c0m0['m0'].mean() + c1m1['m1'].mean())) + 0.5 * cov_c0_inv2 + 0.5 * cov_c1_inv2
float(gamma412)
Out[17]:
1108.7886957742055
In [8]:
#最优化方法——消费计算 使用调整之后的报销比例平衡面板
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['r0adjust'] - data['r1adjust']) * data['m0'] + data['premium2015'] - data['premium2018'])[0, 1]
cov_c1_inv2 = np.cov(data['c1']**(-3) / E_c1_inv2, (data['r0adjust'] - data['r1adjust']) * data['m1'] + data['premium2015'] - data['premium2018'])[0, 1]

# 计算 gamma
data['gamma422'] = abs(data['premium2015'] - data['premium2018']) + abs(0.5 * (data['r0adjust'] - data['r1adjust']) * (E_m0 + E_m1)) + 0.5 * cov_c0_inv2 + 0.5 * cov_c1_inv2
gamma422= data['gamma422'].mean()
print(gamma422)
data

1015.7927795030351
Out[8]:
ID c0 m0 householdID communityID c1 m1 gender age marriage ... r1 r0adjust r1adjust policyintergration2015 policyintergration2018 district GDPgrowthrate urban_nbs gamma322 gamma422
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 704.829975 919.486397
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 704.829975 919.486397
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 704.829975 919.486397
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 704.829975 919.486397
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 704.829975 919.486397
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
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 753.326538 729.140131
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 753.326538 729.140131
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 753.326538 729.140131
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 753.326538 729.140131
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 753.326538 729.140131

3887 rows × 28 columns

health-based：合并数据¶

In [10]:
dataa = pd.merge(c0m0, c1m1full,on="ID", how="inner")
dataa
Out[10]:
ID c0 m0 householdID communityID c1 m1 gender age marriage ... premium2018 r0 r1 r0adjust r1adjust policyintergration2015 policyintergration2018 district GDPgrowthrate urban_nbs
0 64033321002 112.216 60.0 640333210 640333 123.670 0.0 0.0 59.0 1.0 ... 180.0 0.6167 0.700 0.544192 0.660903 0.0 1.0 east 0.118005 Rural
1 64033327002 1029.200 6000.0 640333270 640333 1054.764 21000.0 0.0 62.0 1.0 ... 180.0 0.6167 0.700 0.544192 0.660903 0.0 1.0 east 0.118005 Rural
2 64033325001 1062.400 3000.0 640333250 640333 78.020 100.0 1.0 66.0 1.0 ... 180.0 0.6167 0.700 0.544192 0.660903 0.0 1.0 east 0.118005 Rural
3 64033322001 592.620 2000.0 640333220 640333 859.050 17050.0 1.0 63.0 1.0 ... 180.0 0.6167 0.700 0.544192 0.660903 0.0 1.0 east 0.118005 Rural
4 64033330002 2058.400 2000.0 640333300 640333 4025.500 2000.0 1.0 59.0 1.0 ... 180.0 0.6167 0.700 0.544192 0.660903 0.0 1.0 east 0.118005 Rural
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
3882 89676104001 2315.700 840.0 896761040 896761 891.420 8000.0 1.0 61.0 1.0 ... 220.0 0.6500 0.725 0.624462 0.685108 0.0 1.0 west 0.284050 Rural
3883 89676114002 1935.145 300.0 896761140 896761 49.800 0.0 0.0 56.0 0.0 ... 220.0 0.6500 0.725 0.624462 0.685108 0.0 1.0 west 0.284050 Rural
3884 89676118001 2466.096 1000.0 896761180 896761 661.095 3000.0 0.0 73.0 0.0 ... 220.0 0.6500 0.725 0.624462 0.685108 0.0 1.0 west 0.284050 Rural
3885 89676115001 10721.940 800.0 896761150 896761 11638.260 2000.0 1.0 55.0 1.0 ... 220.0 0.6500 0.725 0.624462 0.685108 0.0 1.0 west 0.284050 Rural
3886 89676124001 268.422 500.0 896761240 896761 313.242 101.0 1.0 69.0 0.0 ... 220.0 0.6500 0.725 0.624462 0.685108 0.0 1.0 west 0.284050 Rural

3887 rows × 26 columns

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

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

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

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

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

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

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  hsf15   R-squared:                       0.000
Model:                            OLS   Adj. R-squared:                  0.000
Method:                 Least Squares   F-statistic:                     1.297
Date:                Mon, 29 Dec 2025   Prob (F-statistic):              0.255
Time:                        17:16:27   Log-Likelihood:                -11963.
No. Observations:                3113   AIC:                         2.393e+04
Df Residuals:                    3111   BIC:                         2.394e+04
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         33.6007      0.216    155.587      0.000      33.177      34.024
m0         -1.706e-05    1.5e-05     -1.139      0.255   -4.64e-05    1.23e-05
==============================================================================
Omnibus:                       41.096   Durbin-Watson:                   1.849
Prob(Omnibus):                  0.000   Jarque-Bera (JB):               32.348
Skew:                          -0.167   Prob(JB):                     9.46e-08
Kurtosis:                       2.629   Cond. No.                     1.54e+04
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.54e+04. This might indicate that there are
strong multicollinearity or other numerical problems.
Out[11]:
np.float64(-1.7062284276591698e-05)
In [12]:
#计算dh/dm 18
import pandas as pd
import statsmodels.api as sm
f1 = pd.merge(c1m1,hsf18,on="ID",how="inner")
f1= f1[['m1','hsf18']].copy()
# 删除包含 NaN 或 inf 的行
f1= f1.replace([np.inf, -np.inf], np.nan).dropna()
# 删除包含 0 的行
f1 = f1[(f1['hsf18']!=0) & (f1['m1']!=0)]
# 自变量（X）和因变量（Y）
X = f1['m1']
Y = f1['hsf18']

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

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

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

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

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

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                  hsf18   R-squared:                       0.000
Model:                            OLS   Adj. R-squared:                  0.000
Method:                 Least Squares   F-statistic:                     1.579
Date:                Mon, 29 Dec 2025   Prob (F-statistic):              0.209
Time:                        17:16:27   Log-Likelihood:                -26371.
No. Observations:                5630   AIC:                         5.275e+04
Df Residuals:                    5628   BIC:                         5.276e+04
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         51.5258      0.359    143.352      0.000      50.821      52.230
m1         -1.359e-05   1.08e-05     -1.257      0.209   -3.48e-05    7.61e-06
==============================================================================
Omnibus:                    14923.644   Durbin-Watson:                   1.748
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              441.401
Skew:                          -0.235   Prob(JB):                     1.42e-96
Kurtosis:                       1.711   Cond. No.                     3.42e+04
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 3.42e+04. This might indicate that there are
strong multicollinearity or other numerical problems.
Out[12]:
np.float64(-1.3590936773055398e-05)
In [13]:
import pandas as pd
import numpy as np

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

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

gamma313=abs(city_gender_age_premium_ratio_district15['premium2015'].mean() - city_gender_age_premium_ratio_district15['premium2018'].mean()) + abs(0.5 * (city_gender_age_premium_ratio_district15['r0'].mean() - city_gender_age_premium_ratio_district15['r1'].mean()) * (c0m0['m0'].mean() + c1m1['m1'].mean())) + 0.5 * cov_c0_inv2 + 0.5 * cov_c1_inv2
float(gamma313)
Out[13]:
785.5178011684857
In [14]:
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['gamma323'] = 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
gamma323 = dataa['gamma323'].mean()
print(gamma323)
dataa

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

3887 rows × 27 columns

In [15]:
import pandas as pd
import numpy as np

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

# 计算协方差
cov_c0_inv2 = np.cov((0.019743 * h_m15)/ (E_c0_inv2 * dataa['r0adjust']), (dataa['r0adjust'] - dataa['r1adjust']) * dataa['m0'] + dataa['premium2015'] - dataa['premium2018'])[0, 1]
cov_c1_inv2 = np.cov((0.019743 * h_m18)/ (E_c1_inv2 * dataa['r1adjust']), (dataa['r0adjust'] - dataa['r1adjust']) * dataa['m1'] + dataa['premium2015'] - dataa['premium2018'])[0, 1]

gamma413=abs(city_gender_age_premium_ratio_district15['premium2015'].mean() - city_gender_age_premium_ratio_district15['premium2018'].mean()) + abs(0.5 * (city_gender_age_premium_ratio_district15['r0adjust'].mean() - city_gender_age_premium_ratio_district15['r1adjust'].mean()) * (c0m0['m0'].mean() + c1m1['m1'].mean())) + 0.5 * cov_c0_inv2 + 0.5 * cov_c1_inv2
float(gamma413)
Out[15]:
824.7701821065729
In [16]:
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['r0adjust']), (dataa['r0adjust'] - dataa['r1adjust']) * dataa['m0'] + dataa['premium2015'] - dataa['premium2018'])[0, 1]
cov_c1_inv2 = np.cov((0.019743 * h_m18)/ (E_c1_inv2 * dataa['r1adjust']), (dataa['r0adjust'] - dataa['r1adjust']) * dataa['m1'] + dataa['premium2015'] - dataa['premium2018'])[0, 1]

# 计算 gamma
dataa['gamma423'] = abs(dataa['premium2015'] - dataa['premium2018']) + abs(0.5 * (dataa['r0adjust'] - dataa['r1adjust']) * (E_m0 + E_m1)) + 0.5 * cov_c0_inv2 + 0.5 * cov_c1_inv2
gamma423 = dataa['gamma423'].mean()
print(gamma423)
dataa

731.7742658354025
Out[16]:
ID c0 m0 householdID communityID c1 m1 gender age marriage ... r1 r0adjust r1adjust policyintergration2015 policyintergration2018 district GDPgrowthrate urban_nbs gamma323 gamma423
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 500.242209 635.467883
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 500.242209 635.467883
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 500.242209 635.467883
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 500.242209 635.467883
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 500.242209 635.467883
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
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 548.738772 445.121617
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 548.738772 445.121617
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 548.738772 445.121617
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 548.738772 445.121617
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 548.738772 445.121617

3887 rows × 28 columns

用完全信息法求解¶

In [12]:
import numpy as np
import pandas as pd

# 参数
sigmamale = 3.0
phi_tilde = 0.019743

# 取各列的均值（忽略缺失）
c0_bar  = pd.to_numeric(c0m0["c0"],    errors="coerce").mean(skipna=True)
c1_bar  = pd.to_numeric(c1m1["c1"],    errors="coerce").mean(skipna=True)
h0_bar  = pd.to_numeric(hsf15["hsf15"], errors="coerce").mean(skipna=True)
h1_bar  = pd.to_numeric(hsf18["hsf18"], errors="coerce").mean(skipna=True)

B_bar = (c0_bar**(1 - sigmamale)) + (1 - sigmamale) * phi_tilde * (h0_bar - h1_bar)
cons1_bar = B_bar**(1 / (1 - sigmamale))

gamma311 = c1_bar - cons1_bar
print(gamma311)

1957.6372848184362
In [13]:
import numpy as np
import pandas as pd
# 参数
sigma = 3.0
phi_tilde = 0.019743

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["B_bar"] = (d3["c0"]**(1 - sigma)) + (1 - sigma) * phi_tilde * (d3["hsf15"] - d3["hsf18"])
d3["gamma221"] = d3["c1"] - d3["B_bar"]**(1 / (1 - sigma))

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

1996.3867642406108
Out[13]:
ID c0 m0 c1 m1 hsf15 hsf18 B_bar gamma221
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