求助lee-carter模型R语言代码

我参照的是——孙佳美的《生命表编制理论与实践》进行经典lee-carter模型死亡率预测。数据如下图 我想用1997-2010年死亡率估计2011-2015年的。年龄组有10组。
代码如下：三种参数估计方法——奇异值分解svd，加权最小二乘法WLS，极大似然估计法likelihood。4个地方报错，求解！
library(demography)
library(forecast)
library(xts)
#读取数据
data = read.csv(file = "F://c15total.csv", header = TRUE)
data
break_year<-2010
T<-break_year-1997+1    #用来建模年份
n<-2015-break_year      #用来预测的年份
D_all<-matrix(rep(0,10*19),10,19)    #年龄组数：10组和年数：14
E_all<-matrix(rep(0,10*19),10,19)
m_all<-matrix(rep(0,10*19),10,19)
lnm_all<-matrix(rep(0,10*19),10,19)
for(t in 1:19){
    for(x in 1:10){
      D_all[x,t]<-mydata[x+10*(t-1),3];
      E_all[x,t]<-mydata[x+10*(t-1),4];
      m_all[x,t]<-D_all[x,t]/E_all[x,t];
      lnm_all[x,t]<-log(m_all[x,t])
    }
}

D<-D_all[,1:T];
E<-E_all[,1:T];
m<-m_all[,1:T];
lnm_all<-lnm_all[,1:T];     #建模数据

m_original<-m_all[,(T+1):19]   #预测数据

age<-mydata[1:10,2];
year<-as.numeric(levels(factor(mydata[,1])))
oldyear<-year[1:T];
newyear<-(T+1):(n+T)

#参数估计——加权最小二乘法WLS
a<-c(rep(0,10));
b<-c(rep(1,10));
k<-c(rep(0,T))      #设置初始值
    #建模
RSS<-function(a,b,k){
  L<-matrix(rep(NA,10*T),10,T)
  for(x in 1:10)
    {for(t in 1:T)
      {L[x,t]<-D[x,t]*((lnm[x,t]-a[x]-b[x]*k[t])^2)
       }
  }
  RSS<-sum(L)
  return(RSS)
}
    #迭代a(x)年龄死亡率的均值
aa<-function(a,b,k){
  for(x in 1:10){
    a[x]<-sum(D[X,]*(lnm[x,]-b[x]*k))/sum(D[x,])
  }
  return(a)
}
    #迭代k(t)时间因子
kk<-function(a,b,k){
  for(t in 1:15){
    k[t]<-sum(D[,t]*k*(lnm[x,]-a))/sum(D[,t]*b^2)
  }
  return(b)
}
    #迭代b(x)年龄因子
bb<-function(a,b,k){
  for(x in 1:10){
    b[x]<-sum(D[x,]*k*(lnm[x,]-a[x]))/sum(D[x,]*k^2)
  }
  return(b)
}
    #迭代过程
j=1
diff_RSS<-RSS(a,b,k)     Error in RSS(a, b, k) : object 'lnm' not found
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
}    Error in abs(diff_RSS) : object 'diff_RSS' not found
    #标化
a1<-a;
b1<-b;
k1<-k
c1<-sum(b1);
c2<-sum(k1)/T
a_wls<-a1+c2*b/c1;
b_wls<-b1/c1;
k_wls<-c1*k1-c2


#参数估计——奇异值分解法SVD
age_svd<-c(1:10)
std.m<-demogdata(m,E,ages = age_svd,years = oldyear,type = "mortality",label = "swden",
                 name = "c15",lambda = 0)
swden<-lca(std.m,years = std.m$year,ages = std.m$age,adjust = "dt",restype = "logrates")
Warning messages:
1: In stats::nlm(fred, guess, ...) : NA/Inf被换成最大的正值
2: In stats::nlm(fred, guess, ...) : NA/Inf被换成最大的正值
3: In stats::nlm(fred, guess, ...) : NA/Inf被换成最大的正值
4: In newroot(FUN, guess, ...) : No root exists. Returning closest
5: In newroot(FUN, guess, ...) : No root exists. Returning closest

a_svd<-swden$ax;b_svd<-swden$bx;k_svd<-swden$kt

#导出参数估计数据a(x),b(x),k(t)
age<-as.data.frame(age)
result1<-as.data.frame(cbind(age,a_svd,a_wls,a_likelihood,b_svd,b_wls,b_likelihood));
result2<-as.data.frame(cbind(year[1:T],k_svd,k_wls,k_likelihood))

write.table(result1,"parameters.csv",append=T,sep="",col.names=c("age","a_svd","a_wls",
            "a_likelihood","b_svd","b_wls","b_likelihood"),row.names=F)
Warning message:
In write.table(result1, "parameters.csv", append = T, sep = "",  :
  给文件加列名

write.table(result2,"parameters.csv",append=T,sep="",col.names=c("year","k_svd","k_wls",
            "k_likelihood"),row.names=F)


Warning message:
In write.table(result1, "parameters.csv", append = T, sep = "",  :
  给文件加列名

#代码结束
结果如下

后面还有用参数估计出来的k(t)根据ARIMA模型预测k_predict(t)，然后再和a(x)、b(x)算出死亡率。我这部分都没有代码，求代码啊求代码。真是要跪了。没懂是什么意思，我用R软件算出的只有40岁以上10个年龄组的a(x)、b(x)值和2011-2015年k_predict值，最后怎么拟合出死亡率的预测值啊？
求大神解答！！！！有偿，非数理专业，医学统计出身

第一个error是L[x,t]<-D[x,t]*((lnm[x,t]-a[x]-b[x]*k[t])^2这句代码里把lnm改成lnm_all,这个错误你在三个函数定义里都出错了，改完之后第二个error也就没有了。

我最近也在研究这个，同求代码！

dy949294155 发表于 2017-11-16 15:52
我最近也在研究这个，同求代码！

你研究到哪一步了？

请问您现在研究到哪一步了？出结果了吗

郭晶晶A 发表于 2018-3-15 20:54
请问您现在研究到哪一步了？出结果了吗

看懂原理后，自己用sas编辑剩余的代码，还ok吧。

天堂的狗坚强 发表于 2017-11-15 15:51
我参照的是——孙佳美的《生命表编制理论与实践》进行经典lee-carter模型死亡率预测。数据如下图我想用1997 ...

请问你现在结果出来了吗？方便分享一下源代码吗？谢谢

您好，请问您有孙佳美老师这本书的电子版吗，我现在买不到这本书

浑源县0 发表于 2019-2-28 22:31
您好，请问您有孙佳美老师这本书的电子版吗，我现在买不到这本书

电子版没有，学校图书馆借的，强烈表白苏大独墅湖图书馆，救我毕业论文。

kyl12 发表于 2019-1-9 23:18
请问你现在结果出来了吗？方便分享一下源代码吗？谢谢

最近写毕业论文，在整理，事比较多，等我提交完盲审论文后会公开代码