r for LC

r for LC

#svd分解法
#导入美国数据
USA.mortality<-read.csv("/Users/liuqing/Downloads/死亡率数据库/美国/美国死亡率1950-2017.csv",header=TRUE)
USA.population<-read.csv("/Users/liuqing/Downloads/死亡率数据库/美国/美国人数.csv",header=TRUE)
break_year<-2012
T<-break_year-1950+1  #T=63,1950-2012 63年用来拟合模型
n<-2017-break_year    #n=5,2012-2017用来预测模型
D_all<-matrix(rep(0,100*68),100,68) #计算的死亡人数
E_all<-matrix(rep(0,100*68),100,68) #总人数
m_all<-matrix(rep(0,100*68),100,68) #死亡率
lnm_all<-matrix(rep(0,100*68),100,68) #ln死亡率
year<-as.numeric(levels(factor(USA.population[,1])))# levels(factor筛选出不重复的项，1950-2017
age<-USA.mortality[1:100,2] #0-99岁

#女性1950-2017年的所有数据
for (t in 1:68){
  for(x in 1:100){
    E_all[x,t]<-USA.population[x+100*(t-1),3];  #100行，68列
    m_all[x,t]<-USA.mortality[x+100*(t-1),3];
    D_all[x,t]<-E_all[x,t]*m_all[x,t];   #计算的死亡人数
    lnm_all[x,t]<-log(m_all[x,t])
  }
}

#用于实验的数据
E<-E_all[,1:T]
m<-m_all[,1:T]
lnm<-lnm_all[,1:T]
D<-D_all[,1:T]

#用于预测的数据
m_original<-m_all[,(T+1):68]  #用于预测的死亡率

model_year<-year[1:T]          #1950-2007
fo_year<-year[(T+1):(n+T)]     #2008-2017
#svd
library(demography)
age_svd<-c(0:99) #0到99岁
std.m<-demogdata(m,E,ages=age_svd,years = model_year,type="mortality",label="USA",name="Female",lambda = 0)
USA<-lca(std.m,years=std.m$year,ages = std.m$age,adjust = "dt",restype = "logrates")
a_svd<-USA$ax
b_svd<-USA$bx
k_svd<-USA$kt

#极大似然法
y<-D_all[,1:63]#实际观测到的死亡人数

#极大似然估计
a<-c(rep(0,100))
b<-c(rep(1,100))
k<-c(rep(0,63))

#1
logL<-function(a,b,k){
  L<-matrix(rep(NA,100*63),100,63)
for(x in 1:100){
  for(t in 1:63){
    L[x,t]<-y[x,t]*(a[x]+b[x]*k[t])-E[x,t]*exp(a[x]+b[x]*k[t])#似然函数的值
    }
 }
  logL<-sum(L)
  return(logL)
}
#2.1
yy<-function(a,b,k){
  y_hat<-matrix(rep(0,100*63),100,63)
  for(t in 1:63){
    for(x in 1:100){
      y_hat[x,t]<-E[x,t]*exp(a[x]+b[x]*k[t])# 估计的死亡人数
    }
  }
  return(y_hat)
}
#2.2
aa<-function(a,b,k){
  y_hat<-yy(a,b,k)
  for(x in 1:100){
    a[x]<-a[x]+sum(y[x,]-y_hat[x,])/(sum(y_hat[x,]))
  }
  return(a)
}
#2.3
kk<-function(a,b,k){
  y_hat<-yy(a,b,k)
  for(t in 1:63){
    k[t]<-k[t]+sum((y[,t]-y_hat[,t])*b)/(sum(y_hat[,t]*(b^2)))
  }
  return(k)
}
#2.4
bb<-function(a,b,k){
  y_hat<-yy(a,b,k)
  for(x in 1:100){
    b[x]<-b[x]+sum((y[x,]-y_hat[x,])*k)/(sum(y_hat[x,]*(k^2)))
  }
  return(b)
}
#3
j=1
diff_L<-logL(a,b,k)
while(abs(diff_L)>10^(-10))
{
  logL0<-logL(a,b,k)
  if(j%%3==1){a<-aa(a,b,k)}
  if(j%%3==2){k<-kk(a,b,k)}
  if(j%%3==0){b<-bb(a,b,k)}
  logL1<-logL(a,b,k)
  diff_L<-logL1-logL0 #对数似然函数新的估计值与它的前一个估计值的偏差足够小
  j=j+1
}
#4标准化
a1<-a;b1<-b;k1<-k
c1<-sum(b1);c2<-c1*sum(k1)/63
a_likelihood<-a1+c2*b1/c1
b_likelihood<-b1/c1
k_likelihood<-c1*k1-c2

plot(age,a_likelihood)
plot(age,b_likelihood)
plot(year,k_likelihood)
#5将估计到参数写入表格
age<-as.data.frame(age)
result1<-as.data.frame(cbind(age,a_likelihood,b_likelihood))
result2<-cbind(year,k_likelihood)
write.table(result1,"/Users/liuqing/Downloads/死亡率数据库/美国/女性/估计参数likelihood.csv",append=T,sep="  ",col.names = c("age","a_likelihood","b_likelihood"),row.names = F)
write.table(result2,"/Users/liuqing/Downloads/死亡率数据库/美国/女性/k_likelihood.csv",append=T,sep=" ",col.names = c("year","k_likelihood"),row.names = F)
#加权最小二乘法
#设置初值
a<-c(rep(0,100))
b<-c(rep(1,100))
k<-c(rep(0,63))
#2.编写residuals
RSS<-function(a,b,k){
  L<-matrix(rep(NA,100*63),100,63)
  for(x in 1:100){
    for(t in 1:63){
      L[x,t]=d[x,t]*((lnm[x,t]-a[x]-b[x]*k[t])^2)
    }
  }
  RSS<-sum(L)
  return(RSS)
}
#2.1 set the rule of a(x)
aa<-function(a,b,k){
  for(x in 1:100){
    a[x]<-sum(d[x,]*(lnm[x,]-b[x]*k))/sum(d[x,])
  }
  return(a)
}
#2.2 set the iteration rule of k(t)
kk<-function(a,b,k){
  for(t in 1:63){
    k[t]<-sum(d[,t]*b*(lnm[,t]-a))/sum(d[,t]*b^2)
  }
  return(k)
}
#2.3 set the iteration rule of b(x)
bb<-function(a,b,k){
  for (x in 1:100){
    b[x]<-sum(d[x,]*k*(lnm[x,]-a[x]))/sum(d[x,]*k^2)
  }
  return(b)
}
# iteration process
j=1
diff_RSS<-RSS(a,b,k)
while(abs(diff_RSS)>10^(-10)){
  RSS0<-RSS(a,b,k)
  if(j%%3==1){a<-aa(a,b,k)}
  if(j%%3==2){k<-kk(a,b,k)}
  if(j%%3==0){b<-bb(a,b,k)}
  RSS1<-RSS(a,b,k)
  diff_RSS<-RSS1-RSS0
  j=j+1
}
#standardize the results
a1<-a
b1<-b
k1<-k
c1<-sum(b1)
c2<-c1*sum(k1)/63
a_wls<- a1+c2*b/c1
b_wls<- b1/c1
k_wls<- c1*k1-c2