```{r}
## MORTALITY RATES DATA PREPARATION ##
ageLimit <- 90   #截止到90岁
ageGroups <- ageLimit + 1
timePeriods <- 29 #一共有29年的数据
tHor <- 10  #预测年份
mrRaw <- read.csv("D:/大学资料/毕业论文/实践/死亡率/mortality_data_male.csv", header=TRUE)                                                
colnames(mrRaw) <- c("Age",1994:2022 )               
mrMa <- mrRaw[-1]/100               #对数进行变化，模拟效果更好
mrMatrix <- as.matrix(mrMa)                       
colnames(mrMatrix) <- 1994:2022                                                                               
rownames(mrMatrix) <- 0:ageLimit        

mode(mrMatrix) = "numeric"
mrMatrix <- replace(mrMatrix, mrMatrix == 0, 0.000000001) #为了避免无穷的出现，进行转换
```


```{r}
## BASIC LEE-CARTER METHOD ##
## Fitting procedure ##
a_x = rowMeans(log(mrMatrix))                                                                               
M = log(mrMatrix) - a_x                                                                                               
svdM <- svd(M)                                                                                                                        #利用SVD方法进行分解
b_x <- -svdM$u[,1]                                                                                                       
k_t <- -svdM$d[1]*svdM$v[,1]                                                                       
Mhat <- matrix(nrow = ageGroups, ncol = timePeriods)                               
for (i in 1:timePeriods){                                                                                       
        Mhat[,i] <- a_x + b_x*k_t[i]   #利用估计的参数对原数据进行拟合
}
colnames(Mhat) <- 1994:2022                                                                                       
rownames(Mhat) <- 0:ageLimit
a_x
b_x
k_t
```


```{r}
## Forecasting procedure ##
dhat <- (k_t[timePeriods] - k_t[1])/(timePeriods - 1)                                        # Define drift parameter
total = matrix(nrow = timePeriods - 1, ncol = 1)                                       
for (i in 1:timePeriods-1) {                                                                                       
        total[i] <- (k_t[i+1] - k_t[i] - dhat)^2
}
see = sqrt(1/(timePeriods-2)*sum(total))                                                               
kForc_t <- matrix(nrow = tHor, ncol = 1)                                                                # 估计k_t向量

for (i in 1:tHor) {                                                                                                                #
        kForc_t[i] = k_t[timePeriods] + dhat*i + sqrt(i)*rnorm(1,0,see^2)        # Random walk
}
Mtilde <- matrix(nrow = ageGroups, ncol = tHor)                                                        # Initiate M-tilde matrix
for (i in 1:tHor) {                                                                                                                # 估计未来数据
        Mtilde[,i] <- a_x + b_x*kForc_t[i]
}
colnames(Mtilde) <- 2023:2032                                                                               
rownames(Mtilde) <- 0:ageLimit
Mstar <- cbind(Mhat,Mtilde)                                                                                                # 进行合并
```