预期寿命

求预期寿命

#LC
#导入生命表
USA.lifetableF<-read.csv("/Users/liuqing/Downloads/死亡率数据库/美国/美国女性生命表1933-2017.csv",header=TRUE)
lx<-matrix(rep(0,100*85),100,85)   #真实的生存人数1933-2017
ex<-matrix(rep(0,100*85),100,85)   #真实的余命
qx<-matrix(rep(0,100*85),100,85)   #实际观测的死亡人数


year<-1933:2017#85年
age<-0:99      #100个
cohort<-1851:2012 #100+63-1=162

#女性1950-2017年的所有数据
for (t in 1:85){
  for(x in 1:100){
    lx[x,t]<-USA.lifetableF[x+100*(t-1),6];   
    qx[x,t]<-USA.lifetableF[x+100*(t-1),4];    
    ex[x,t]<-USA.lifetableF[x+100*(t-1),10] 
  }
}
plot(year,ex[66,],col="blue",cex=0.5,lty=2,ann=TRUE,xlab = "year",ylab = "e65",xaxt="n")
axis(1,seq(1933,2017,5))

#预测值
#预测k_t(用泊松分布的似然法)
#k_likelihood是1950:2012（63年的）
library(nlme)
library(forecast)
library(xts)
library(ggplot2)
arima1<-auto.arima(k_likelihood,trace=T)       #Best model: ARIMA(0,1,0) with drift 
fit1<-arima(k_likelihood,order=c(0,1,0),seasonal=list(order=c(1,1,1),period=12))#
f.k<-forecast(fit1,h=14,level=c(99.5))#置信区间为99.5%

Point Forecast Lo 99.5 Hi 99.5
64 -41.74018 -46.75340 -36.72697
65 -43.52865 -50.61841 -36.43889
66 -45.34533 -54.02847 -36.66219
67 -46.42219 -56.44861 -36.39576
68 -46.76886 -57.97874 -35.55897
69 -48.87187 -61.15169 -36.59206
70 -50.49467 -63.75839 -37.23096
71 -51.05032 -65.22982 -36.87082
72 -53.27584 -68.31547 -38.23621
73 -54.46015 -70.31331 -38.60699
74 -55.30361 -71.91155 -38.69567
75 -56.60523 -73.93511 -39.27536
76 -58.62139 -77.10572 -40.13706
77 -60.54980 -80.12062 -40.97899

k_predict<-c(-41.74018,-43.52865,-45.34533,-46.42219,-46.76886,-48.87187,-50.49467,-51.05032,-53.27584,-54.46015,-55.30361,-56.60523,-58.62139,-60.54980)
k_new<-c(k_likelihood,k_predict) #1950:2026年的kt
m_yuce<-matrix(rep(0,100*77),100,77)  #1950:2026（77年的）
for( x in 1:100){
  for(t in 1:77){
    m_yuce[x,t]<-exp(a_likelihood[x]+b_likelihood[x]*k_new[t])
  }
}
miu<-matrix(rep(0,100*77),100,77)
q<-matrix(rep(0,100*77),100,77)
miu<-m_yuce
for( x in 1:100){
  for(t in 1:77){
    q[x,t]<-1-exp(-miu[x,t])
  }
}
l<-lx[,18:80]#1950-2012
#还要预测 2013-2026
l_yuce<-matrix(rep(0,100*14),100,14)
for( x in 1:100){
  for(t in 1:14){
    l_yuce[x,t]<-E[x,t]*(1-exp(a_likelihood[x]+b_likelihood[x]*k_predict[t]))
  }
}
#将预测值跟实际值放在一个矩阵里
l_new<-cbind(l,l_yuce)#100*77

#计算余命(65岁时的)
e_shouming<-matrix(rep(0,100*77),100,77)
  #for(t in 1:77){
    for(i in 0:34){
      e_shouming<-(sum(l_new[66+i,1]*(1-(1/2)*q[66+i,1])))/l[66,1]
    }
#}
e_shouming