import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
import matplotlib.gridspec as gridspec
%matplotlib inline
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
from mpl_toolkits.mplot3d import Axes3D
import seaborn as sns
import statsmodels.formula.api as smf
import statsmodels.api as sm
import copy
import warnings
warnings.filterwarnings("ignore")
from sklearn.linear_model import LinearRegression
from scipy.stats import gaussian_kde
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error
%%time
def do_rawdata(df0): #函数注释：从长护险试点原始数据中提取状态转移数据
    df = df0.copy()
    max_year = df.disable_date.max().year
    max_month = df.disable_date.max().month
    data_l = []
    for i_ in  df.id.unique():
        d0 = df[df.id==i_].sort_values('age').reset_index(drop=True)
        d1 = clean_p_data(d0,max_year,max_month)
        d2  = get_chgdata(d1)
        data_l.append(d2)
    chg_data =pd.concat(data_l,axis=0)
    return chg_data[chg_data.delta>=0].reset_index(drop=True)
def clean_p_data(d,max_year,max_month):#函数注释：判断样本是否在数据截止日，并添加最后时刻健康状态
    res = pd.DataFrame()
    age = d.age.iloc[0]
    id_ = d['id'].iloc[0]
    l_level = d.level.to_list()
    l_dt = pd.to_datetime(d.disable_date)
    t0 = l_dt.iloc[0]
    l_dt2 = l_dt.apply(lambda x:(x-t0).days).to_list()
    if d.dtime.isna().values[0]:
        l_level.append(l_level[-1])
        t = (pd.to_datetime('{}-{}-01'.format(max_year,max_month+1)) - l_dt.iloc[0]).days
        l_dt2.append(t)
    else:
        l_level.append('4')#与正文一致，状态4表示死亡
        t = (pd.to_datetime(d.dtime.values[0]) - l_dt.iloc[0]).days
        if t<0:
            l_dt2.append(30)
        else:
            l_dt2.append(t)
        
    res['state'] = l_level
    res['id'] = id_
    res['age'] = age
    res['delta'] = l_dt2
    res['date'] = l_dt.sort_values().iloc[0]
    return res[['id','age','delta','state','date']]
def get_chgdata(d_):#函数注释：提取月度层面连续变化的健康状态数据,具体过程见附录V
    n = len(d_)
    l_s = d_.state.values
    l_dt = d_.delta.apply(lambda x:int(x/30)+1)
    ag = d_.age.iloc[0] *12
    id_ = d_.id.values[0]
    l_ag1 = []
    l_dt2 = []
    l_s1 = []
    l_s2 = []
    n = len(l_s)
    for i in range(n-1):
        s1 = l_s[i]
        s2 = l_s[i+1]
        ag1 = ag+l_dt[i]
        dt = l_dt[i+1]-l_dt[i]
        l_ag1.append(ag1)
        l_dt2.append(dt)
        l_s1.append(s1)
        l_s2.append(s2)
    res = pd.DataFrame()
    res['age1'] = l_ag1
    res['delta'] =l_dt2
    res['state1'] = l_s1
    res['state2'] = l_s2
    res['date'] = d_.date.values[0]
    res['id'] = d_.id.values[0]
    return res

df = pd.read_excel('某副省级城市长护险试点数据.xlsx')
df = df[df.disable_date.apply(lambda x:len(str(x))>4)]#删除建档时间为空
df.disable_date = pd.to_datetime(df.disable_date)#建档时间提取为时间对象
df['year'] = df.disable_date.apply(lambda x:x.year)
df = df[df.level.notna()]#保留所有失能存在失能等级且建档的样本
df.level = df.level.apply(lambda x:x[-1])#提取失能等级的级数
chgdata =  do_rawdata(df)



def get_px(t1,d,s1):#函数注释：计算单时刻状态转移概率
    dic = {'1':0,'2':0,'3':0,'4':0}
    d0 = d.copy()
    d0['age2'] = d0.age1+d0.delta
    d1 = d0[(d0.age1<=t1) & (d0.age2>t1+1)]
    d2 = d0[(d0.age1<=t1) & (d0.age2>t1) & (d0.age2<=(t1+1))]
    n = len(d1)+len(d2)
    if n>0:
        dic[s1] = len(d1)+dic[s1]
        sr = d2.state2.value_counts()
        for i in range(1,5):
            if str(i) in sr.index:
                dic[str(i)] = dic[str(i)]+ sr[str(i)]
            dic[str(i)] = dic[str(i)]/n

        return list(dic.values()),n
    else:
        return [np.nan,np.nan,np.nan,np.nan],0
    
def get_px_res(df,s1):#函数注释：整理状态转移概率数据
    d = df[df.state1 ==s1] 
    l0  = [] 
    l1 = []
    l2 = []
    l3 = []
    l4 = []
    n_l = []
    for  i in np.arange(600,1261,1):
        l,n_ = get_px(i,d,s1)
        l0.append(i)
        l1.append(l[0])
        l2.append(l[1])
        l3.append(l[2])
        l4.append(l[3])
        n_l.append(n_)
    r = pd.DataFrame()
    r['age'] = l0
    r['px1'] = l1
    r['px2'] = l2
    r['px3'] = l3
    r['px4'] = l4
    r['n'] = n_l
    return r

def func_px(t,p,w):#函数注释：通过加权最小二乘方法对状态转移概率进行拟合
    bol = p.notna()
    w0 = w[bol]
    w0 = w0.apply(lambda x:x/w0.sum())
    X = pd.DataFrame()
    X['t'] = t
    X['t2'] = t.apply(lambda x:x**2)
    X['t3'] = t.apply(lambda x:x**3) 
    X['t4'] = t.apply(lambda x:x**4) 
    X = sm.add_constant(X)
    model = LinearRegression()
    model.fit(X[bol],p[bol],sample_weight=w0)
    
    return [model.predict(X),model]


def get_matrix(df):#函数注释：提取状态转移矩阵
    r1 = get_px_res(df,'1')
    r2 = get_px_res(df,'2')
    r3 = get_px_res(df,'3')
    t = r1.age
    w1,w2,w3 = r1.n, r2.n, r3.n
    r = pd.DataFrame()
    r['age'] = r1.age
    r['p11'] ,r['p12'],r['p13'],r['p14'] = func_px(t,r1.px1,w1)[0],func_px(t,r1.px2,w1)[0],func_px(t,r1.px3,w1)[0],func_px(t,r1.px4,w1)[0]
    r['p21'] ,r['p22'],r['p23'],r['p24'] = func_px(t,r2.px1,w2)[0],func_px(t,r2.px2,w2)[0],func_px(t,r2.px3,w2)[0],func_px(t,r2.px4,w2)[0]
    r['p31'] ,r['p32'],r['p33'],r['p34'] = func_px(t,r3.px1,w3)[0],func_px(t,r3.px2,w3)[0],func_px(t,r3.px3,w3)[0],func_px(t,r3.px4,w3)[0]
    r['p41'] ,r['p42'],r['p43'],r['p44'] = 0,0,0,1
    
    r = r.set_index('age')
    r = correct(r).reset_index()
    return r
def correct(m):#函数注释：对状态转移矩阵进行修正
    l = ['p'+str(i)+str(j) for i in range(1,4) for j in range(1,5)]
    res = m.copy()
    for l0 in l:
        s = res[l0].copy()
        if s.min()<0:
            s = s-s.min()
        for i in range(1,len(s)+1):
            if s.iloc[-i]<0:
                try:
                    s.iloc[-i] = s.iloc[-i+1]*s.iloc[-i+1]/s.iloc[-i+2]
                except:
                    continue
        res[l0] = s
    s1 = res.iloc[:,:4].sum(axis=1)
    s2 = res.iloc[:,4:8].sum(axis=1)
    s3 = res.iloc[:,8:12].sum(axis=1)
    for l0 in l:
        s = res[l0].copy()
        if l0[1]=='1':
            s = s/s1
        elif l0[1] =='2':
            s = s/s2
        elif l0[1]=='3':
            s = s/s3
        res[l0] = s
    return res


fig = plt.figure(figsize=(30,24),dpi=800)
#========================================
r1 = get_px_res(chgdata,'1').fillna(0)
r2 = get_px_res(chgdata,'2').fillna(0)
r3 = get_px_res(chgdata,'3').fillna(0)
t = r1.age
w1,w2,w3 = r1.n, r2.n, r3.n
cmap = plt.cm.viridis
r = pd.DataFrame()
r['age'] = r1.age
r['p11'] ,r['p12'],r['p13'],r['p14'] = func_px(t,r1.px1,w1)[0],func_px(t,r1.px2,w1)[0],func_px(t,r1.px3,w1)[0],func_px(t,r1.px4,w1)[0]
r['p21'] ,r['p22'],r['p23'],r['p24'] = func_px(t,r2.px1,w2)[0],func_px(t,r2.px2,w2)[0],func_px(t,r2.px3,w2)[0],func_px(t,r2.px4,w2)[0]
r['p31'] ,r['p32'],r['p33'],r['p34'] = func_px(t,r3.px1,w3)[0],func_px(t,r3.px2,w3)[0],func_px(t,r3.px3,w3)[0],func_px(t,r3.px4,w3)[0]


def get_text(y,y_pred,w0):#函数注释：批注拟合效果
    r2 = r2_score(y, y_pred, sample_weight=w0)
    mse = mean_squared_error(y, y_pred, sample_weight=w0)
    mae = mean_absolute_error(y, y_pred, sample_weight=w0)
    rmse = np.sqrt(mse)

    # 将统计量整理成 Series 对象
    stats_series = pd.Series({
        '可决系数(r2_score)': r2,
        '均方误差 (MSE)': mse,
        '平均绝对误差 (MAE)': mae,
        '根均方误差 (RMSE)': rmse
    })
    text = ""
    for index, value in stats_series.items():
        text += f"{index}: {value:.4f}\n"
    return text
def plot_reg(w0,y_pred,y,c_value,range_,ax,v):#函数注释：拟合过程的可视化函数
    plt.bar(x=idx[::12],height=c_value,width=11,color=cmap(norm(c_value)),label='拟合权重(N(1,t))',alpha=0.6,zorder=1)
    plt.xticks([])
    plt.ylim(0,range_[0])
    plt.ylabel('样本量')
    handles2, labels2 = ax.get_legend_handles_labels()
    ax2 = ax.twinx()
    sns.scatterplot(x=idx,y=y,color=sns.color_palette()[1],marker = 'o',s=8,label='Prb({}\u2192{})'.format(v[1],v[2]),alpha=1,zorder=0)
    plt.ylim(range_[1],range_[2])
    plt.xlabel('年龄t')
    plt.ylabel('转移概率')
    ax2.set_xticks(idx[::12]) 
    ax2.set_xticklabels([int(i/12) for i in r.age[::12]],fontsize=7)
    sns.lineplot(x=idx,y=y_pred,label='拟合结果')
    handles, labels = ax2.get_legend_handles_labels()
    plt.xticks(rotation=45)
    text = get_text(y,y_pred,w0)
    plt.xlim(120,600)
    plt.text(130, range_[1], text, fontsize=20)
    plt.legend(handles+handles2,labels+labels2,fontsize=10,loc='upper right')
##——————————————————————————————————————————状态转移概率P11的拟合结果——————————————————————————————————————
norm = plt.Normalize(min(w1), max(w1))
idx = [i for i in range(len(r))]
w0 = w1/w1.sum()
v = 'p11'
y_pred = r[v]
y = r1.px1
c_value = r1.n[::12]
range_ = [1800,0.9,1.01]
ax = fig.add_subplot(4,3,1)
plot_reg(w0,y_pred,y,c_value,range_,ax,v)
##——————————————————————————————————————————状态转移概率P12的拟合结果——————————————————————————————————————
ax = fig.add_subplot(4,3,2)
w0 = w1/w1.sum()
v = 'p12'
y_pred = r[v]
y = r1.px2
c_value = r1.n[::12]
range_ = [1800,-0.02,0.02]
plot_reg(w0,y_pred,y,c_value,range_,ax,v)
##——————————————————————————————————————————状态转移概率P13的拟合结果——————————————————————————————————————
ax = fig.add_subplot(4,3,3)
w0 = w1/w1.sum()
v = 'p13'
y_pred = r[v]
y = r1.px3
c_value = r1.n[::12]
range_ = [1800,-0.02,0.02]
plot_reg(w0,y_pred,y,c_value,range_,ax,v)
##——————————————————————————————————————————状态转移概率P14的拟合结果——————————————————————————————————————
ax = fig.add_subplot(4,3,4)
w0 = w1/w1.sum()
v = 'p14'
y_pred = r[v]
y = r1.px4
c_value = r1.n[::12]

range_ = [1800,0,0.06]
plot_reg(w0,y_pred,y,c_value,range_,ax,v)

##——————————————————————————————————————————状态转移概率P21的拟合结果——————————————————————————————————————
w0 = w2/w2.sum()
v = 'p21'
y_pred = r[v]
y = r2.px1
c_value = r2.n[::12]
ax = fig.add_subplot(4,3,5)
range_ = [1500,-0.01,0.01]
plot_reg(w0,y_pred,y,c_value,range_,ax,v)
##——————————————————————————————————————————状态转移概率P22的拟合结果——————————————————————————————————————
w0 = w2/w2.sum()
v = 'p22'
y_pred = r[v]
y = r2.px2
c_value = r2.n[::12]
ax = fig.add_subplot(4,3,6)
range_ = [1500,0.9,1]
plot_reg(w0,y_pred,y,c_value,range_,ax,v)
##——————————————————————————————————————————状态转移概率P23的拟合结果——————————————————————————————————————
w0 = w2/w2.sum()
v = 'p23'
y_pred = r[v]
y = r2.px3
c_value = r2.n[::12]
ax = fig.add_subplot(4,3,7)
range_ = [1500,-0.01,0.01]
plot_reg(w0,y_pred,y,c_value,range_,ax,v)
##——————————————————————————————————————————状态转移概率P24的拟合结果——————————————————————————————————————
w0 = w2/w2.sum()
v = 'p24'
y_pred = r[v]
y = r2.px4
c_value = r2.n[::12]
ax = fig.add_subplot(4,3,8)
range_ = [1500,0.0,0.06]
plot_reg(w0,y_pred,y,c_value,range_,ax,v)
##——————————————————————————————————————————状态转移概率P33的拟合结果——————————————————————————————————————
w0 = w3/(w3.sum())
v = 'p33'
y_pred = r[v]
y = r3.px3
c_value = r3.n[::12]
ax = fig.add_subplot(4,3,9)
range_ = [20,0.6,1.1]
plot_reg(w0,y_pred,y,c_value,range_,ax,v)
##——————————————————————————————————————————状态转移概率P34的拟合结果——————————————————————————————————————
w0 = w3/(w3.sum())
v = 'p34'
y_pred = r[v]
y = r3.px4
c_value = r3.n[::12]
ax = fig.add_subplot(4,3,10)
range_ = [20,-0.1,0.5]
plot_reg(w0,y_pred,y,c_value,range_,ax,v)
plt.show()

def f_m(age,r):#获得单个矩阵
    m = r[r.age==age].values[0][1:].reshape(4,4)
    return m 
def get_int_matrix(r):#获得整数矩阵
    res = pd.DataFrame()
    l_ = np.arange(600,1261,12)
    for i in range(len(l_)-1):
        a0 = l_[i]
        an = l_[i+1]
        m_ = np.eye(4)
        for j in range(a0,an+1):
            m_ = np.dot(m_,f_m(j,r))
        res[a0] = [p[0] for p in list(m_.reshape(16,1))]
    res = res.T.reset_index()
    res.columns = r.columns
    res.age = res.age/12
    return res

def get_all_matrix(r):#获得全部矩阵
    res = []
    l_ = np.arange(600,1261)
    for i in range(len(l_)-1):
        # print(i)
        a1 = l_[i]
        m_ = np.eye(4)
        for a2 in range(a1+1,1261):
            m_ = np.dot(m_,f_m(a2-1,r))
            res0 = pd.DataFrame()
            res0[a2] = [p[0] for p in list(m_.reshape(16,1))]
            res0 = res0.T.reset_index()
            res0.columns = r.columns
            res0['delta'] = a2-a1
            res0['age'] = a1
            res.append(res0)
    return pd.concat(res)
r = get_matrix(chgdata)
matrix = get_all_matrix(r)

def plot_swt_xy(r0,a0,p_name): #函数注释：提取状态转移概率
    r = r0[(r0.age==a0) | (r0.age==a0+1)]
    x = np.arange(a0,a0+2)
    n = r.delta.max()
    y = np.arange(1,n)
    r = r[r.delta<n].sort_values(['delta','age'])
    x,y = np.meshgrid(x, y)
    n_x = x.shape[1]
    n_y = y.shape[0]
    lx = x.reshape(1,n_x*n_y)[0]
    ly = y.reshape(1,n_x*n_y)[0]
    z1 = np.array(r[p_name].values.reshape(n_y,n_x))
    return x,y,z1

def swt(ax,x,y,z,xl='',yl='',zl='',el=20,az=150,cmap ='binary'): #函数注释：用于转移概率模型的可视化
    ax.plot_surface(x,y,z,rstride=1, cstride=1,alpha=1, cmap=plt.get_cmap(cmap))
    ax.set_ylabel(yl)
    ax.set_xlabel(xl)
    ax.set_zlabel(zl)
fig = plt.figure(figsize=(24,15),dpi=600)
#——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,1,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p11')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200)  
plt.title('P11(t,'+chr(916)+'t)',y=0.8)
#——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,2,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p12')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200)  
plt.title('P12(t,'+chr(916)+'t)',y=0.8)
# #——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,3,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p13')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200) 
plt.title('P13(t,'+chr(916)+'t)',y=0.8)
# #——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,4,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p14')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200)  
plt.title('P14(t,'+chr(916)+'t)',y=0.8)
# #——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,5,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p21')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200) 
plt.title('P21(t,'+chr(916)+'t)',y=0.8)
# #——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,6,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p22')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200) 
plt.title('P22(t,'+chr(916)+'t)',y=0.8)
# #——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,7,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p23')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200) 
plt.title('P23(t,'+chr(916)+'t)',y=0.8)
# #——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,8,projection="3d")

for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p24')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200) 
plt.title('P24(t,'+chr(916)+'t)',y=0.8)
# #——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,9,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p33')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200)  
plt.title('P33(t,'+chr(916)+'t)',y=0.8)
# #——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,10,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p34')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200)  
plt.title('P34(t,'+chr(916)+'t)',y=0.8)
plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.05, hspace=-0.05)

def plot_swt_xy(r0,a0,p_name): #函数注释：提取状态转移概率
    r = r0[(r0.age==a0) | (r0.age==a0+1)]
    x = np.arange(a0,a0+2)
    n = r.delta.max()
    y = np.arange(1,n)
    r = r[r.delta<n].sort_values(['delta','age'])
    x,y = np.meshgrid(x, y)
    n_x = x.shape[1]
    n_y = y.shape[0]
    lx = x.reshape(1,n_x*n_y)[0]
    ly = y.reshape(1,n_x*n_y)[0]
    z1 = np.array(r[p_name].values.reshape(n_y,n_x))
    return x,y,z1

def swt(ax,x,y,z,xl='',yl='',zl='',el=20,az=150,cmap ='binary'): #函数注释：用于转移概率模型的可视化
    ax.plot_surface(x,y,z,rstride=1, cstride=1,alpha=1, cmap=plt.get_cmap(cmap))
    ax.set_ylabel(yl)
    ax.set_xlabel(xl)
    ax.set_zlabel(zl)
fig = plt.figure(figsize=(24,15),dpi=600)
#——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,1,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p11')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200)  
plt.title('P11(t,'+chr(916)+'t)',y=0.8)
#——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,2,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p12')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200)  
plt.title('P12(t,'+chr(916)+'t)',y=0.8)
# #——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,3,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p13')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200) 
plt.title('P13(t,'+chr(916)+'t)',y=0.8)
# #——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,4,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p14')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200)  
plt.title('P14(t,'+chr(916)+'t)',y=0.8)
# #——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,5,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p21')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200) 
plt.title('P21(t,'+chr(916)+'t)',y=0.8)
# #——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,6,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p22')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200) 
plt.title('P22(t,'+chr(916)+'t)',y=0.8)
# #——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,7,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p23')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200) 
plt.title('P23(t,'+chr(916)+'t)',y=0.8)
# #——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,8,projection="3d")

for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p24')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200) 
plt.title('P24(t,'+chr(916)+'t)',y=0.8)
# #——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,9,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p33')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200)  
plt.title('P33(t,'+chr(916)+'t)',y=0.8)
# #——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,4,10,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix,a,'p34')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p')
ax.view_init(elev=10,azim=200)  
plt.title('P34(t,'+chr(916)+'t)',y=0.8)
plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.05, hspace=-0.05)

def get_chgyear(df0):#函数注释：将总状态转移数据分解为不同年度状态数据
    df = df0.copy()
    df['year'] = df.date.apply(lambda x:x.year)
    df['month'] = df.date.apply(lambda x:x.month)+1
    df['i'] = df.month+df.delta
    df_1 = df[df.i<=13].reset_index(drop=True)
    df_2 = df[(df.i>13)&(df.i<=25) ].reset_index(drop=True)
    df_3 = df[(df.i>25)&(df.i<=37) ].reset_index(drop=True)
    df_4 = df[(df.i>37)].reset_index(drop=True)
    res = pd.concat([decompose(df_2,2),decompose(df_3,3),decompose(df_4,4)])
    df_1['year'] = df_1.date.apply(lambda x:x.year)
    df_1 = df_1[res.columns]
    return pd.concat([res,df_1])
def decompose(df0,t):#函数注释：将总状态转移数据分解为不同年度状态数据
    df = df0.copy()
    df1,df2,df3,df4 = pd.DataFrame(),pd.DataFrame(),pd.DataFrame(),pd.DataFrame()
    df1,df2 = pd.DataFrame(),pd.DataFrame()
    df1['age1'] = df.age1
    df1['delta'] = df.month.apply(lambda x:13-x)
    df1['state1'],df1['state2'] = df.state1,df.state1
    df1['year'] = df.year
    df1['id'] = df.id
    
    df2['age1'] = df1.delta+df.age1
    df2['delta'] = df.delta-df1.delta
    df2['state1'],df2['state2'] =df.state1,df.state2
    df2['year'] =df.year+1
    df2['id'] = df.id
    if t==2:
        df1['date'],df2['date'] = df.date,df.date
        return pd.concat([df1,df2],axis=0)
    elif t==3:
        df3['age1'] = df1.delta+12+df.age1
        df3['delta'] = df.delta-df1.delta-12
        df3['state1'],df3['state2'] =df.state1,df.state2
        df3['year'] =df.year+2
        df3['id'] = df.id
        df1['date'],df2['date'] ,df3['date']= df.date,df.date,df.date
        return  pd.concat([df1,df2,df3],axis=0)
    elif t==4:
        df3['age1'] = df1.delta+df.age1+12
        df3['delta'] = 12
        df3['state1'],df3['state2'] =df.state1,df.state1
        df3['year'] =df.year+2
        df3['id'] = df.id

        df4['age1'] = df1.delta+df.age1+24
        df4['delta'] = df.delta-df1.delta-24
        df4['state1'],df4['state2'] =df.state1,df.state2
        df4['year'] =df.year+3
        df4['id'] = df.id

        df1['date'],df2['date'] ,df3['date'],df4['date']= df.date,df.date,df.date,df.date
        return pd.concat([df1,df2,df3,df4],axis=0)

chgdata_year = get_chgyear(chgdata)
chgdata_year_17 = get_matrix(chgdata_year[chgdata_year.year==2017])
matrix_17 = get_all_matrix(chgdata_year_17)
chgdata_year_18 = get_matrix(chgdata_year[chgdata_year.year==2018])
matrix_18 = get_all_matrix(chgdata_year_18) 
chgdata_year_19 = get_matrix(chgdata_year[chgdata_year.year==2019])
matrix_19 = get_all_matrix(chgdata_year_19)
chgdata_year_20 = get_matrix(chgdata_year[chgdata_year.year==2020])
matrix_20 = get_all_matrix(chgdata_year_20)

def get_lp(r,age,s,delta): #函数注释：提取单个状态转移矩阵，计算思路与正文一致
    r0 = r[(r.age==age)&(r.delta==delta)]
    l_p = r0[['p{}1'.format(s),'p{}2'.format(s),'p{}3'.format(s),'p{}4'.format(s)]].values[0]
    return l_p

def get_vec_q(r,age,s,delta): #函数注释：计算状态向量，计算思路与正文一致
    for i in range(delta):
        if i==0:
            res = get_lp(r,age,s,delta=1)
            l_s = res
        else:
            res2= res[0]*get_lp(r,age+i,s=1,delta=1)+res[1]*get_lp(r,age+i,s=2,delta=1)+res[2]*get_lp(r,age+i,s=3,delta=1)+res[3]*get_lp(r,age+i,s=4,delta=1)
            res=res2
            l_s += res2
    return l_s
def rest_data(r): #函数注释：计算不同状态转移矩阵下的失能余寿
    idx = [i for i in range(720,1260,12)]
    l_1 = [get_vec_q(r,age=i,s=1,delta=1260-i) for i in idx]
    l_2 = [get_vec_q(r,age=i,s=2,delta=1260-i) for i in idx]

    res = pd.DataFrame()
    res['age'] = idx
    res.age = res.age/12
    res['r1'] = [i[:3].sum()/12 for i in l_1]
    res['r2'] = [i[:3].sum()/12 for i in l_2]
    res['vec1'] = l_1
    res['vec2'] = l_2
    return res

rest_age = rest_data(matrix)     #计算平均余寿
rest_age17 = rest_data(matrix_17)#计算2017平均余寿
rest_age18 = rest_data(matrix_18)#计算2018平均余寿
rest_age19 = rest_data(matrix_19)#计算2019平均余寿
rest_age20 = rest_data(matrix_20)#计算2020平均余寿

#================================计算第七次全国人口普查数据中的平均余寿============================
d0 =pd.read_excel('6-4 全国分年龄性别的死亡人口状况2019.11.1-2020.10.31.xls',header=3)
d_city =pd.read_excel('6-4a 全国分年龄性别的死亡人口状况2019.11.1-2020.10.31城市.xls',header=3)
d_rural =pd.read_excel('6-4c 全国分年龄性别的死亡人口状况2019.11.1-2020.10.31乡村.xls',header=3)
def res_age(data,age):
    d = data[0]
    d_rate_100=data[1]
    
    d0= d[d.age>=age]
    s = d0.a_rate
    s = pd.concat([s,pd.Series([1-d_rate_100]*15)])#默认115岁以上死亡率为1
    l = []
    for n in range(0,len(s)):
        s2 = s.iloc[0:n+1]
        d_rate = 1-s2.prod()
        l.append(d_rate)
    res = pd.Series(l).diff()
    res[0] = l[0]
    res = sum([(i+1)*res.iloc[i] for i in range(len(s))])
    return res
def do(d):
    d = d[['年  龄','平均人口','死亡人口','死亡率']]
    d.columns = ['age','num','d_num','d_rate']
    d_rate_100 = d.iloc[-1,:]['d_rate']/1000
    d =d[d.age.apply(lambda x:str(x) in [str(i) for i in range(100)])]
    d = d.astype('f')
    d['d_rate'] = d.d_num/d.num
    d['a_rate'] = 1-d.d_rate
    return [d,d_rate_100]
res = pd.DataFrame()
res['age'] = [i for i in range(60,101)]
res['r'] = [res_age(do(d0),i) for i in range(60,101)] #计算全国人口的平均余寿
res['r_city'] = [res_age(do(d_city),i) for i in range(60,101)]#计算全国城市人口的平均余寿
res['r_rural'] = [res_age(do(d_city),i) for i in range(60,101)]#计算全国乡村人口的平均余寿

#=====================================与长护险试点数据中的平均余寿进行比较=======================
fig = plt.figure(figsize=(20,6),dpi=600)
ax = fig.add_subplot(1,2,1)

plt.plot(rest_age.age,rest_age.r1,marker='o',color=sns.color_palette()[0],label='期望余寿：重度失能一级')
plt.plot(rest_age.age,rest_age.r2,marker='o',color=sns.color_palette()[1],label='期望余寿：重度失能二级')
plt.xlim(59,100)
plt.legend(fontsize=13.6)
plt.grid()
ax = fig.add_subplot(1,2,2)

plt.plot(rest_age.age,rest_age.r1,marker='o',color=sns.color_palette()[0],label='期望余寿：长护险试点数据（重度失能1级）')
plt.plot(rest_age.age,rest_age.r2,marker='o',color=sns.color_palette()[1],label='期望余寿：长护险试点数据（重度失能2级）')
plt.plot(res.age,res.r,marker='o',color=sns.color_palette()[2],label='期望余寿：第七次全国人口普查数据（全部人口）')
plt.plot(res.age,res.r_city,marker='o',color=sns.color_palette()[4],label='期望余寿：第七次全国人口普查数据（城市人口）')
plt.plot(res.age,res.r_rural,marker='o',color=sns.color_palette()[5],label='期望余寿：第七次全国人口普查数据（农村人口）')


plt.legend(fontsize=13.6)
plt.grid()
plt.xlim(59,100)
plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.05, hspace=0.2)

%%time
fig = plt.figure(figsize=(24,30),dpi=600)
ft=10
#——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(6,4,1,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix_17[matrix_17.delta<360],a,'p14')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p',cmap='Reds' )
ax.view_init(elev=12,azim=170) 
plt.title('2017年重度1级死亡率',y=0.9,fontsize=ft)
#——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(6,4,2,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix_18[matrix_18.delta<360],a,'p14')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p',cmap='Wistia')
ax.view_init(elev=12,azim=170) 
plt.title('2018年重度1级死亡率',y=0.9,fontsize=ft)
#——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(6,4,5,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix_19[matrix_19.delta<360],a,'p14')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p',cmap='BuGn')
ax.view_init(elev=12,azim=170) 
plt.title('2019年重度1级死亡率',y=0.9,fontsize=ft)
#——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(6,4,6,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix_20[matrix_20.delta<360],a,'p14')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p',cmap='PuBu')
ax.view_init(elev=12,azim=170) 
plt.title('2020年重度1级死亡率',y=0.9,fontsize=ft)
#——————————————————————————————————————————————————————————————————————————————
# #——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,2,3,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix_17[matrix_17.delta<360],a,'p14')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p',cmap='Reds' )
    x,y,z = plot_swt_xy(matrix_18[matrix_18.delta<360],a,'p14')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p',cmap='Wistia')
    x,y,z = plot_swt_xy(matrix_19[matrix_19.delta<360],a,'p14')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p',cmap='BuGn')
    x,y,z = plot_swt_xy(matrix_20[matrix_20.delta<360],a,'p14')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p',cmap='PuBu')
ax.view_init(elev=12,azim=170) 
plt.title('各年度重度1级死亡率对比',y=0.9,fontsize=ft+10)
#————————————————————————————————————————————————————————————————————————————————————————————
#——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(6,4,3,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix_17[matrix_17.delta<360],a,'p24')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p',cmap='Reds' )
ax.view_init(elev=12,azim=170) 
plt.title('2017年重度2级死亡率',y=0.9,fontsize=ft)
#——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(6,4,4,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix_18[matrix_18.delta<360],a,'p24')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p',cmap='Wistia')
ax.view_init(elev=12,azim=170)
plt.title('2018年重度2级死亡率',y=0.9,fontsize=ft)
#——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(6,4,7,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix_19[matrix_19.delta<360],a,'p24')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p',cmap='BuGn')
ax.view_init(elev=12,azim=170) 
plt.title('2019年重度2级死亡率',y=0.9,fontsize=ft)
#——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(6,4,8,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix_20[matrix_20.delta<360],a,'p24')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p',cmap='PuBu')
ax.view_init(elev=12,azim=170) 
plt.title('2020年重度2级死亡率',y=0.9,fontsize=ft)
#——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(3,2,4,projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix_17[matrix_17.delta<360],a,'p24')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p',cmap='Reds' )
    x,y,z = plot_swt_xy(matrix_18[matrix_18.delta<360],a,'p24')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p',cmap='Wistia')
    x,y,z = plot_swt_xy(matrix_19[matrix_19.delta<360],a,'p24')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p',cmap='BuGn')
    x,y,z = plot_swt_xy(matrix_20[matrix_20.delta<360],a,'p24')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p',cmap='PuBu')
ax.view_init(elev=12,azim=170)
plt.title('各年度重度2级死亡率对比',y=0.9,fontsize=ft+10)

ax = fig.add_subplot(3,2,5)

plt.plot(rest_age17.age,rest_age17.r1,marker='o',color=sns.color_palette()[3],label='2017年期望余寿：重度失能1级')
plt.plot(rest_age18.age,rest_age18.r1,marker='o',color=sns.color_palette()[1],label='2018年期望余寿：重度失能1级')
plt.plot(rest_age19.age,rest_age19.r1,marker='o',color=sns.color_palette()[2],label='2019年期望余寿：重度失能1级')
plt.plot(rest_age20.age,rest_age20.r1,marker='o',color=sns.color_palette()[0],label='2020年期望余寿：重度失能1级')
plt.xlim(60,100)
plt.legend(fontsize=13.6)
plt.grid()
ax = fig.add_subplot(3,2,6)
plt.plot(rest_age17.age,rest_age17.r2,marker='o',color=sns.color_palette()[3],label='2017年期望余寿：重度失能2级')
plt.plot(rest_age18.age,rest_age18.r2,marker='o',color=sns.color_palette()[1],label='2018年期望余寿：重度失能2级')
plt.plot(rest_age19.age,rest_age19.r2,marker='o',color=sns.color_palette()[2],label='2019年期望余寿：重度失能2级')
plt.plot(rest_age20.age,rest_age20.r2,marker='o',color=sns.color_palette()[0],label='2020年期望余寿：重度失能2级')
plt.xlim(60,100)
plt.legend(fontsize=13.6)
plt.grid()

plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.05, hspace=-0.05)

sz = pd.read_excel('某副省级城市长护险试点数据.xlsx')
sz = sz[sz.disable_date.apply(lambda x:len(str(x))>4)]#删除建档时间为空
sz.disable_date = pd.to_datetime(sz.disable_date)#建档时间提取为时间对象
sz['year'] = sz.disable_date.apply(lambda x:x.year)
sz  =sz[sz.level.isna()]#所有非失能建档样本即失智样本
sz.level='1'

def get_matrix(df):#函数注释：由于失智样本的健康状态数为2，对get_matrix函数作针对性调整
    r1 = get_px_res(df,'1')
    t = r1.age
    w1= r1.n
    r = pd.DataFrame()
    r['age'] = r1.age
    r['p11'] ,r['p12'],r['p13'],r['p14'] = func_px(t,r1.px1,w1)[0],0,0,func_px(t,r1.px4,w1)[0]
    r['p21'] ,r['p22'],r['p23'],r['p24'] = 0,1,0,0
    r['p31'] ,r['p32'],r['p33'],r['p34'] = 0,0,1,0
    r['p41'] ,r['p42'],r['p43'],r['p44'] = 0,0,0,1
    
    r = r.set_index('age')
    r = correct(r).reset_index()
    return r

chgdata_sz = do_rawdata(sz)
chgdata_sz_year = get_chgyear(chgdata_sz)
chgdata_sz_year_19 = get_matrix(chgdata_sz_year[chgdata_sz_year.year==2019])
chgdata_sz_year_20 = get_matrix(chgdata_sz_year[chgdata_sz_year.year==2020])

matrix_sz_19 = get_all_matrix(chgdata_sz_year_19) 
matrix_sz_20 = get_all_matrix(chgdata_sz_year_20) 

rest_sz_19 = rest_data(matrix_sz_19)
rest_sz_20 = rest_data(matrix_sz_20)

%%time
fig = plt.figure(figsize=(30,10),dpi=600)
spec = gridspec.GridSpec(1, 4, figure=fig)
ft=10
# #——————————————————————————————————————————————————————————————————————————————
ax = fig.add_subplot(spec[0,:1])
plt.plot(rest_sz_19.age,rest_sz_19.r1,marker='o',color=sns.color_palette()[2],label='2019年期望余寿：失智')
plt.plot(rest_sz_20.age,rest_sz_20.r1,marker='o',color=sns.color_palette()[0],label='2020年期望余寿：失智')
plt.legend(fontsize=13.6)
plt.grid()
plt.xlim(65,100)
ax = fig.add_subplot(spec[0,1:],projection="3d")
for a in np.arange(600,1200):
    x,y,z = plot_swt_xy(matrix_sz_19[matrix_sz_19.delta<360],a,'p14')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p',cmap='BuGn')
    x,y,z = plot_swt_xy(matrix_sz_20[matrix_sz_20.delta<360],a,'p14')
    swt(ax,x/12,y/12,z,'年龄t','转移间隔'+chr(916)+'t','转移概率p',cmap='PuBu')
ax.view_init(elev=12,azim=170) 
plt.title('各年度失智群体的死亡率对比',y=0.9,fontsize=ft)

plt.subplots_adjust(left=0.5, right=1, bottom=0.5, top=1, wspace=-0.4, hspace=0.2)

df = pd.read_excel('某副省级城市长护险试点数据.xlsx')
df = df[df.disable_date.apply(lambda x:len(str(x))>4)]#删除建档时间为空
df.disable_date = pd.to_datetime(df.disable_date)#建档时间提取为时间对象
df['year'] = df.disable_date.apply(lambda x:x.year)
df = df[df.level.notna()]#保留所有失能存在失能等级且建档的样本
df.level = df.level.apply(lambda x:x[-1])#提取失能等级的级数

d = df[df.level.isin(['1','2'])].drop_duplicates('id')
d = d[(d.age>=60)&(d.age<100)]
rng = pd.date_range('2017-9','2027-12',freq = 'm')

def cal_pay_abs(r,age,s,year0,month0,year,month):
    delta = (year-year0)*12+month-month0  
    r0 = r[(r.age==age*12)&(r.delta==delta)]
    l_p = r0[['p{}1'.format(s),'p{}2'.format(s),'p{}3'.format(s),'p{}4'.format(s)]].values
    if l_p.shape[0]>0:
        l_p2 = l_p[0]
        return l_p2[0]*0.6+l_p2[1]*0.8+l_p2[2]
    return 0
def cal_pay_abstime(r,d,year,month):
    d0=d.copy()
    d0['year0']  =d0.disable_date.apply(lambda x:x.year)
    d0['month0']  =d0.disable_date.apply(lambda x:x.month)
    d1 = d0[['age','level','year0','month0']].drop_duplicates()
    d1['pay'] = d1.apply(lambda x:cal_pay_abs(r=r,age=x.age,s=x.level,year0=x.year0,month0=x.month0,year=year,month=month),axis=1)
    d2 = pd.merge(d1,d0,on=['age','level','year0','month0'],how='right')
    s = d2['pay']
    return s.sum()
def cal_age_abs(r,age,s,year0,month0,year,month):
    delta = (year-year0)*12+month-month0
    r0 = r[(r.age==age*12)&(r.delta==delta)]
    l_p = r0[['p{}1'.format(s),'p{}2'.format(s),'p{}3'.format(s),'p{}4'.format(s)]].values
    if l_p.shape[0]>0:
        l_p2 = l_p[0]
        return l_p2[0]+l_p2[1]+l_p2
    return 0
def cal_age_abstime(r,d,year,month):
    d0=d.copy()
    d0['year0']  =d0.disable_date.apply(lambda x:x.year)
    d0['month0']  =d0.disable_date.apply(lambda x:x.month)
    d1 = d0[['age','level','year0','month0']].drop_duplicates()
    d1['pay'] = d1.apply(lambda x:cal_age_abs(r=r,age=x.age,s=x.level,year0=x.year0,month0=x.month0,year=year,month=month),axis=1)
    d2 = pd.merge(d1,d0,on=['age','level','year0','month0'],how='right')
    s = d2['pay']
    return s.sum()

res = pd.DataFrame()
res['dt'] = rng
def cal(r,d,x):
    year = x.year
    month = x.month
    return cal_pay_abstime(r,d,year,month)
res['pay17'] = res.dt.apply(lambda x:cal(matrix_17,d,x))
res['pay18'] = res.dt.apply(lambda x:cal(matrix_18,d,x))
res['pay19'] = res.dt.apply(lambda x:cal(matrix_19,d,x))
res['pay20'] = res.dt.apply(lambda x:cal(matrix_20,d,x))
res.to_pickle('pay绝对时间基金支出.pkl')

def cal_pay(r,age,s,delta):
    r0 = r[(r.age==age*12)&(r.delta==delta)]
    l_p = r0[['p{}1'.format(s),'p{}2'.format(s),'p{}3'.format(s),'p{}4'.format(s)]].values
    if l_p.shape[0]>0:
        l_p2 = l_p[0]
        return l_p2[0]*0.6+l_p2[1]*0.8+l_p2[2]
    return 0
def cal_pay_time(r,d,delta):
    d0=d.copy()
    d1 = d0[['age','level']].drop_duplicates() 
    d1['pay'] = d1.apply(lambda x:cal_pay(r=r,age=x.age,s=x.level,delta=delta),axis=1)
    d2 = pd.merge(d1,d0,on=['age','level'],how='right')
    s = d2['pay']
    return s.sum()
def cal_age(r,age,s,delta):
    r0 = r[(r.age==age*12)&(r.delta==delta)]
    l_p = r0[['p{}1'.format(s),'p{}2'.format(s),'p{}3'.format(s),'p{}4'.format(s)]].values
    if l_p.shape[0]>0:
        l_p2 = l_p[0]
        return l_p2[0]+l_p2[1]+l_p2[2]
    return 0
def cal_age_time(r,d,delta):
    d0=d.copy()

    d1 = d0[['age','level']].drop_duplicates()
    
    d1['pay'] = d1.apply(lambda x:cal_age(r=r,age=x.age,s=x.level,delta=delta),axis=1)
    d2 = pd.merge(d1,d0,on=['age','level'],how='right')
    s = d2['pay']
    return s.sum()

res2 = pd.DataFrame()
res2['delta'] = [i+1 for i in range(len(rng))]
res2['pay17'] = res2['delta'].apply(lambda x: cal_pay_time(matrix_17,d,x))
res2['pay18'] = res2['delta'].apply(lambda x: cal_pay_time(matrix_18,d,x))
res2['pay19'] = res2['delta'].apply(lambda x: cal_pay_time(matrix_19,d,x))
res2['pay20'] = res2['delta'].apply(lambda x: cal_pay_time(matrix_20,d,x))

res2['age17'] = res2['delta'].apply(lambda x: cal_age_time(matrix_17,d,x))
res2['age18'] = res2['delta'].apply(lambda x: cal_age_time(matrix_18,d,x))
res2['age19'] = res2['delta'].apply(lambda x: cal_age_time(matrix_19,d,x))
res2['age20'] = res2['delta'].apply(lambda x: cal_age_time(matrix_20,d,x))
res2.to_pickle('pay相对时间基金支出age期望余寿.pkl')

fig = plt.figure(figsize=(20,12),dpi=600)
ax = fig.add_subplot(2,2,1)
idx = [i for i in range(len(res))]

plt.plot(idx,res.pay17*1796/10000,marker='o',color=sns.color_palette()[3],label='基金支付水平：基于2017年的平均健康状态',markersize=3)
plt.plot(idx,res.pay18*1796/10000,marker='o',color=sns.color_palette()[1],label='基金支付水平：基于2018年的平均健康状态',markersize=3)
plt.plot(idx,res.pay19*1796/10000,marker='o',color=sns.color_palette()[2],label='基金支付水平：基于2019年的平均健康状态',markersize=3)
plt.plot(idx,res.pay20*1796/10000,marker='o',color=sns.color_palette()[0],label='基金支付水平：基于2020年的平均健康状态',markersize=3)

plt.fill_between((0,38),(2000,2000),color='grey',alpha=0.5)
ax.set_xticks([i for i in idx][::6])  # 设置刻度标签位置
ax.set_xticklabels(['{}-{}'.format(i.year,i.month) for i in res.dt][::6],fontsize=10,rotation=30)
plt.legend(fontsize=10)
plt.ylabel('长护险基金的月度支出水平（万元）',fontsize=15)

ax = fig.add_subplot(2,2,2)
idx = [i for i in range(len(res))]

plt.plot(idx,res.pay17.cumsum()*1796/10000,marker='o',color=sns.color_palette()[3],label='基金累计支付水平：基于2017年的平均健康状态',markersize=3)
plt.plot(idx,res.pay18.cumsum()*1796/10000,marker='o',color=sns.color_palette()[1],label='基金累计支付水平：基于2018年的平均健康状态',markersize=3)
plt.plot(idx,res.pay19.cumsum()*1796/10000,marker='o',color=sns.color_palette()[2],label='基金累计支付水平：基于2019年的平均健康状态',markersize=3)
plt.plot(idx,res.pay20.cumsum()*1796/10000,marker='o',color=sns.color_palette()[0],label='基金累计支付水平：基于2020年的平均健康状态',markersize=3)

plt.fill_between((0,38),(130000,130000),color='grey',alpha=0.5)
ax.set_xticks([i for i in idx][::6])  # 设置刻度标签位置
ax.set_xticklabels(['{}-{}'.format(i.year,i.month) for i in res.dt][::6],fontsize=10,rotation=30)
plt.legend(fontsize=10)
plt.ylabel('长护险基金的月度支出水平（万元）',fontsize=15)

ax = fig.add_subplot(2,2,3)
idx = [i for i in range(len(res2))]


plt.plot(idx,(res2.age17)/len(d),marker='o',color=sns.color_palette()[3],label='平均存活比例：基于2017年的平均健康状态',markersize=3)
plt.plot(idx,(res2.age18)/len(d),marker='o',color=sns.color_palette()[1],label='平均存活比例：基于2018年的平均健康状态',markersize=3)
plt.plot(idx,(res2.age19)/len(d),marker='o',color=sns.color_palette()[2],label='平均存活比例：基于2019年的平均健康状态',markersize=3)
plt.plot(idx,(res2.age20)/len(d),marker='o',color=sns.color_palette()[0],label='平均存活比例：基于2020年的平均健康状态',markersize=3)

ax.set_xticks([i for i in idx][::6])  # 设置刻度标签位置
ax.set_xticklabels(['试点第{}个月'.format(i) for i in res2.delta][::6],fontsize=10,rotation=30)
plt.grid()
plt.ylabel('平均存活比例',fontsize=15)

ax = fig.add_subplot(2,2,4)
v = [res2.age17.sum()/len(d),res2.age18.sum()/len(d),res2.age19.sum()/len(d),res2.age20.sum()/len(d)]
plt.bar([2017,2018,2019,2020],v,width=0.8,edgecolor='black'
        ,color=[sns.color_palette()[3],sns.color_palette()[1],sns.color_palette()[2],sns.color_palette()[0]])

ax.set_xticks([2017,2018,2019,2020])  # 设置刻度标签位置
ax.set_xticklabels(['{}年'.format(i) for i in [2017,2018,2019,2020]],fontsize=10)
plt.ylabel('平均寿命(月)',fontsize=15)