#
#
#  This file is to demonstrate the graduation of mortality rates 
#    in Taiwan 2004 
#  --> 1st and 2nd columns: Male populations and number of deaths 
#      3rd and 4th columns: Female populations and number of deaths 
#
data=matrix(scan("e:\\Taiwan2004.txt"),ncol=4,byrow=T) 
age=c(0:20)*5
#
#  (1) MWA graduation:
#  
mwa=function(data,coef) {
 n=nrow(data)
  m=(length(coef)-1) 
   ux=data[,2]
    vx=rep(0,n) 
     for (i in (m+1):(n-m)) { 
      for (j in (-m:m)) { 
       vx[i]=vx[i]+coef[abs(j)+1]*ux[i+j]
       }
      }
     k=length(vx) 
   vx[1:m]=NA 
  vx[(k-m+1):k]=NA    
 return(vx)
} 
#
#
coef1=c(.4857,.3429,-.0857)
coef2=c(.5594,.2937,-.0734)
#coef1=c(.2554,.2338,.1688,.0606,-.0909)
#coef2=c(.3311,.2666,.1185,-.0099,-.0407)
#
#
# We shall try the male first!
ux=data[,2]/data[,1]
 v1=mwa(cbind(age,ux),coef1)
  v2=mwa(cbind(age,ux),coef2)
matplot(cbind(age,age,age),log(cbind(ux,v1,v2)),type="l",col=1,lty=c(3,1,4),xlab="Age")
#
#
#
#  Note: Since the differences are small, we will use the data in Table 5.3 
#        (You could try data in Table 5.2 as well.)
#
data=matrix(scan("e:\\g53.txt"),ncol=4,byrow=T) 
 age=data[,1]
#
ux=data[,2] 
 v1=mwa(cbind(age,ux),coef1)
  v2=mwa(cbind(age,ux),coef2)
matplot(cbind(age,age,age),cbind(ux,v1,v2),type="l",col=c(2,1,1),lty=c(3,1,4),xlab="Age")
#########################################################################################
#
# (2) Whitakker graduation:
#
wh=function(data,W,z,h=0.1)
{
  n=nrow(data)	# the length of the data
  age=data[,1]
  ux=data[,2]
  K=matrix(rep(0,n*(n-z)),n-z,n)
  for (j  in  1: (n-z))  {
    for ( i  in  0:z)  {
     K[j,i+j]=(-1)^(z-i)*gamma(z+1)/(gamma(z-i+1)*gamma(i+1))  
      } 
     }
  vx=solve((W+h*crossprod(K,K)))%*%W%*%ux
  return(vx)
}  
#
#  Note: The weight here is 1 for all age. 
data=matrix(scan("e:\\g53.txt"),ncol=4,byrow=T) 
x1=data[,1]
 ux=data[,2]
  W=diag(nrow(data))
   v1=wh(data,W,2,h=0.1)
    v2=wh(data,W,2,h=1)
   v3=wh(data,W,2,h=10)
  v4=wh(data,W,2,h=100)
matplot(cbind(x1,x1,x1,x1,x1),cbind(ux,v1,v2,v3,v4),type="l",col=c(2,1,1,3,3),lty=c(3,1,4),xlab="Age")
#
#
#  Here is the program for Taiwan 2004 
data=matrix(scan("g:\\Taiwan2004.txt"),ncol=4,byrow=T) 
 age=c(0:20)*5
  ux=data[,2]/data[,1] 
   W=diag(data[,1]) 
data=cbind(data[,1],ux)
   v1=wh(data,W,2,h=0.1)
    v2=wh(data,W,2,h=1)
   v3=wh(data,W,2,h=10)
  v4=wh(data,W,2,h=100)
matplot(cbind(age,age,age,age,age),cbind(ux,v1,v2,v3,v4),type="l",col=c(2,1,1,3,3),lty=c(3,1,4),xlab="Age")