R for LC+cohort

R for LC+cohort

#导入美国数据
USA.mortality<-read.csv("/Users/liuqing/Downloads/死亡率数据库/美国/美国死亡率1950-2017.csv",header=TRUE)
USA.population<-read.csv("/Users/liuqing/Downloads/死亡率数据库/美国/美国人数1950-2017.csv",header=TRUE)
USA.deathpopulation<-read.csv("/Users/liuqing/Downloads/死亡率数据库/美国/美国死亡人数1950-2017.csv",header=TRUE)
USA.gdp<-read.csv("/Users/liuqing/Downloads/死亡率数据库/美国/女性/美国人均gdp.csv",head=TRUE)#1960-2018
m_all<-matrix(rep(0,100*68),100,68)   #死亡率
lnm_all<-matrix(rep(0,100*68),100,68) #ln死亡率
D_all<-matrix(rep(0,100*68),100,68)   #计算的死亡人数
E_all<-matrix(rep(0,100*68),100,68)   #总人数
d_all<-matrix(rep(0,100*68),100,68)   #实际观测的死亡人数
#女性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]<-USA.deathpopulation[x+100*(t-1),3] #实际观测的死亡人数
    lnm_all[x,t]<-log(m_all[x,t]);
  }
}
#用于实验的数据
E<-E_all[,1:63] #1950-2012
m<-m_all[,1:63] #1950-2012
lnm<-lnm_all[,1:63]
m_original<-lnm_all[,64:68]#用于预测的死亡率
y<-D_all[,1:63]#实际观测到的死亡人数
#添加age-period-cohort
a<-c(rep(0,100))
b0<-c(rep(1,100))
b1<-c(rep(1,100))
k<-c(rep(0,63))
l<-c(rep(0,162))
#1.1 估计的死亡率
mm<-function(a,b0,b1,k,l){
  m_hat<-matrix(rep(0,100*63),100,63)
  for(t in 1950:2012){
    for(x in 0:99){
      for(z in 1851:2012){
        if(z==t-x){
          m_hat[x+1,t-1949]<-exp(a[x+1]+b0[x+1]*k[t-1949]+b1[x+1]*l[t-x-1850])
        }
      }
    }
  }
  return(m_hat)
}
#1.2 估计的死亡人数
yy<-function(a,b0,b1,k,l){
  y_hat<-matrix(rep(0,100*63),100,63)
  E<-E_all[,1:63]
  for(x in 1:100){
    for(t in 1:63){
        y_hat[x,t]<-E[x,t]*m_hat[x,t]      #估计的死亡人数
      }
  }
  return(y_hat)
}


#1.3 似然函数的值
RSS<-function(a,b0,b1,k,l){
  L<-matrix(rep(NA,100*63),100,63)
  y_hat<-yy(a,b0,b1,k,l)
  m_hat<-mm(a,b0,b1,k,l)
  for(x in 1:100){
    for(t in 1:63){
          L[x,t]<-2*(y[x,t]*log(y[x,t]/y_hat[x,t])-(y[x,t]-y_hat[x,t]))
    }
    }
  suml<-sum(L)
  return(suml)
}

#2.1 alpha
aa<-function(a,b0,b1,k,l){
  y_hat<-yy(a,b0,b1,k,l)
  for(x in 1:100){
    a[x]<-a[x]+sum(y[x,]-y_hat[x,])/(sum(y_hat[x,]))
  }
  return(a)
}

#2.2 kt
kk<-function(a,b0,b1,k,l){
  y_hat<-yy(a,b0,b1,k,l)
  for(t in 1:63){
    k[t]<-k[t]+sum((y[,t]-y_hat[,t])*b1)/(sum(y_hat[,t]*(b1^2)))
  }
  return(k)
}

#2.3 beta(1)
bb1<-function(a,b0,b1,k,l){
  y_hat<-yy(a,b0,b1,k,l)
  for(x in 1:100){
    b1[x]<-b1[x]+sum((y[x,]-y_hat[x,])*k)/(sum(y_hat[x,]*(k^2)))
  }
  return(b1)
}

#2.4 cohort
ll<-function(a,b0,b1,k,l){
  y_hat<-yy(a,b0,b1,k,l) 
  for(x in 0:99){           #最后对x循环
    for(z in 1851:2012){    #再对z循环
      for(t in 1950:2012){  #先对t循环
         if(t-x==z){
          l[z-1850]<-l[z-1850]+sum((y[x+1,t-1949]-y_hat[x+1,t-1949])*b0[x+1])/(sum(y_hat[x+1,t-1949]*(b0[x+1]^2)))
        }
      }
    }
  }
  return(l)
}

#2.5
bb0<-function(a,b0,b1,k,l){
  y_hat<-yy(a,b0,b1,k,l)
  for( x in 0:99){
    for(t in 1950:2012){
      for(z in 1851:2012){
        if(z==t-x){
          b0[x+1]<-b0[x+1]+sum((y[x+1,t-1949]-y_hat[x+1,t-1949])*l[t-x-1850])/sum(y_hat[x+1,t-1949]*(l[t-x-1850])^2)
        }
      }
  return(b0)
    }
  }
}
#3
j=1
diff_L<-RSS(a,b0,b1,k,l)
while(abs(diff_L)>10^(-10))
{
  logL0<-RSS(a,b0,b1,k,l)
  if(j%%5==1){a<-aa(a,b0,b1,k,l)}
  if(j%%5==2){k<-kk(a,b0,b1,k,l)}
  if(j%%5==3){b0<-bb0(a,b0,b1,k,l)}
  if(j%%5==4){b1<-bb1(a,b0,b1,k,l)}
  if(j%%5==0){l<-ll(a,b0,b1,k,l)}
  logL1<-RSS(a,b0,b1,k,l)
  diff_L<-logL1-logL0
  j=j+1
}