age+GDP参数估计

#调用美国女性数据
library(demography)
USA <- hmd.mx("USA", "alkaid139@163.com", "19970618", "USA")
#
USA2<-lifetable(USA,years=1950:2012,age=0:99,type="period")#生命表
USA3<-hmd.pop("USA", "alkaid139@163.com", "19970618", "USA")#总人数
USA.death<-read.csv("//HN/HPCUserFile$/liuqing/User/Desktop/美国死亡人数1950-2017.csv",header=TRUE)#1960-2018
#death 
d_all<-matrix(rep(0,100*68),100,68) 
#女性1950-2017年的所有数据
for (t in 1:68){
  for(x in 1:100){
    d_all[x,t]<-USA.death[x+100*(t-1),3] #实际观测的死亡人数
  }
}
d<-d_all[,11:63]#1960-2012
E_female<-((USA3$pop)[1])#
E_female$female
E<-(E_female$female)[1:100,28:80]

USA.gdp<-read.csv("//HN/HPCUserFile$/liuqing/User/Desktop/美国人均gdp的副本.csv",header=TRUE)#1960-2018
summary(USA.gdp)
g_all<-log(USA.gdp$USA人均.GDP.2010年不变价美元.)
g<-g_all[1:53]#1960-2012
#普通最小二乘估计
lnm<-log((USA2$mx)[,11:63]) #100*53 1960-2012
for(x in 1:100){
  a_hat[x]<-(sum(lnm[x,]))/53
}
RSS<-function(b,r,k){
  for(x in 1:100){
    for(t in 1:53){
      L<-(lnm[x,t]-a_hat[x]-b[x]*k[t]-r[x]*g[t])^2
    }
  }
  RSS<-sum(L)
  return(RSS)
}
bb<-function(b,r,k){
  for(x in 1:100){
    b[x]<-sum((lnm[x,]-a_hat[x]-r[x]*g)*k)/(sum(k^2))
  }
  return(b)
}
rr<-function(b,r,k){
  for(x in 1:100){
    r[x]<-sum((lnm[x,]-a_hat[x]-b[x]*k)*g)/(sum(g^2))
  }
  return(r)
}
kk<-function(b,r,k){
  for(t in 1:53){
    k[t]<-sum((lnm[,t]-a_hat-r*g[t])*b)/(sum(b^2))
  }
  return(k)
}
b<-c(rep(1,100))
k<-c(rep(0,53))
r<-c(rep(0,100))
diff_L<-RSS(b,r,k)
diff_L
while(abs(diff_L)>10^(-10))
{
  logL0<-RSS(b,r,k)
  k<-kk(b,r,k)
  r<-rr(b,r,k)
  b<-bb(b,r,k)
  logL1<-RSS(b,r,k)
  diff_L<-logL1-logL0 #对数似然函数新的估计值与它的前一个估计值的偏差足够小
}
d<-sum(b)
e<-cov(k,g)/var(g)
c<-(-sum(k-e*g)/53)
a_hat<-a_hat-b*c
b_hat<-b/d
k_hat<- d*(k-e*g+c)
r_hat<- r+e*b

plot(a_hat)
plot(b_hat)
plot(r_hat)
plot(k_hat)
age<-as.data.frame(0:99)
result1<-as.data.frame(cbind(age,a_hat,b_hat,r_hat))
result2<-cbind(1960:2012,k_hat)
write.table(result1,"//HN/HPCUserFile$/liuqing/User/Desktop/abr11估计参数.csv",append=T,sep="  ",col.names = c("age","a_hat","b_hat","r_hat"),row.names = F)
write.table(result2,"//HN/HPCUserFile$/liuqing/User/Desktop/k估计参数.csv",append=T,sep="  ",col.names = c("year","k_hat"),row.names = F)

#加权估计
RSS1<-function(a,b,r,k){
  for(x in 1:100){
    for(t in 1:53){
      L<-(d[x,t]*(lnm[x,t]-a[x]-b[x]*k[t]-r[x]*g[t]))^2
    }
  }
  RSS1<-sum(L)
  return(RSS1)
}
aa1<-function(a,b,r,k){
  for(x in 1:100){
    a[x]<-sum(d[x,]*(lnm[x,]-b[x]*k-r[x]*g))/sum(d[x,])
  }
  return(a)
}
bb1<-function(a,b,r,k){
  for(x in 1:100){
    b[x]<-sum(d[x,]*k*(lnm[x,]-a[x]-r[x]*g))/(sum(d[x,]*(k^2)))
  }
  return(b)
}
rr1<-function(a,b,r,k){
  for(x in 1:100){
    r[x]<-sum(d[x,]*(lnm[x,]-a[x]-b[x]*k)*g)/(sum(d[x,]*(g^2)))
  }
  return(r)
}
kk1<-function(a,b,r,k){
  for(t in 1:53){
    k[t]<-sum(d[,t]*(lnm[,t]-a-r*g[t])*b)/(sum(d[,t]*(b^2)))
  }
  return(k)
}
a<-c(rep(1,100))
b<-c(rep(1,100))
k<-c(rep(0,53))
r<-c(rep(0,100))
diff_L<-RSS1(a,b,r,k)
diff_L
while(abs(diff_L)>10^(-10))
{
  logL0<-RSS1(a,b,r,k)
  k<-kk1(a,b,r,k)
  r<-rr1(a,b,r,k)
  b<-bb1(a,b,r,k)
  logL1<-RSS1(a,b,r,k)
  diff_L<-logL1-logL0 #对数似然函数新的估计值与它的前一个估计值的偏差足够小
}
d<-sum(b)
e<-cov(k,g)/var(g)
c<-(-sum(k-e*g)/53)


a_hat<-a-b*c
b_hat<-b/d
k_hat<- d*(k-e*g+c)
r_hat<- r+e*b

plot(b_hat)
plot(r_hat)
plot(k_hat)
age<-as.data.frame(0:99)
result1<-as.data.frame(cbind(age,a_hat,b_hat,r_hat))
result2<-cbind(1960:2012,k_hat)
write.table(result1,"//HN/HPCUserFile$/liuqing/User/Desktop/abr1估计参数.csv",append=T,sep="  ",col.names = c("age","a_hat","b_hat","r_hat"),row.names = F)
write.table(result2,"//HN/HPCUserFile$/liuqing/User/Desktop/k1估计参数.csv",append=T,sep="  ",col.names = c("year","k_hat"),row.names = F)

#热度图
canshu<-read.csv("/Users/liuqing/Downloads/远程桌面文件/abr估计参数.csv",header=TRUE)
kkkk<-read.csv("/Users/liuqing/Downloads/远程桌面文件/k估计参数.csv",header=TRUE)
a<-canshu$a_hat
b<-canshu$b_hat
r<-canshu$r_hat
k<-kkkk$k_hat

a1<-canshu$a_hat.weight
b1<-canshu$b_hat.Weight
r1<-canshu$r_hat.weight
k1<-kkkk$k_hat.weight


e<-matrix(rep(0,100*53),100,53)
lnm1<-lnm[,11:63]
for(x in 1:100){
  for(t in 1:53){
    e[x,t]<-lnm1[x,t]-a[x]-b[x]*k[t]-r[x]*g[t]
  }
}

e1<-matrix(rep(0,100*53),100,53)
lnm1<-lnm[,11:63]
for(x in 1:100){
  for(t in 1:53){
    e1[x,t]<-sqrt(D[x,t])*(lnm1[x,t]-a1[x]-b1[x]*k1[t]-r1[x]*g[t])
  }
}

#画热度图
library(pheatmap)
library(RColorBrewer)
colnames(e)=c(1960:2012)
rownames(e)=c(0:99)
pheatmap(e,
         color=brewer.pal(11,"PRGn"),
         cellwidth = 3,#每个单元格的宽度
         cellheight = 3,#每个单元格的高度
         cluster_rows = FALSE,
         cluster_cols = FALSE,
         labels_row=c(0:99),
         fontsize = 3,#热图中字体显示的大小
         #display_numbers = TRUE,
         number_color = "black"
)
#画热度图
colnames(e1)=c(1960:2012)
rownames(e1)=c(0:99)
pheatmap(e1,
         color=brewer.pal(11,"PRGn"),
         cellwidth = 3,#每个单元格的宽度
         cellheight = 3,#每个单元格的高度
         cluster_rows = FALSE,
         cluster_cols = FALSE,
         labels_row=c(0:99),
         fontsize = 2,#热图中字体显示的大小
         #display_numbers = TRUE,
         number_color = "black"
)