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