nesau — Mar 12, 2014, 10:45 AM
## This program matches the output in the back of DHW
A = 0.00022
B = 2.7e-06
c = 1.124
x = 20
w = 131
i = 0.05
l20 = 100000
formatVector <- function(nums, k) {
for (j in 1:length(nums)) {
num = as.numeric(nums[j])
nums[j] = format(round(num,k), nsmall=k)
}
return(as.numeric(nums))
}
formatDataFrame <- function(table, numDecimals) {
for (j in (1:length(table))) {
numEntries = w-x
currentcol = table[0:numEntries,][j][0:numEntries,]
currentcol = formatVector(currentcol, numDecimals)
currentcol = data.frame(currentcol)
table[0:numEntries,][j] = currentcol
}
return(table)
}
display <- function(message, num) {
cat(paste(message, num, sep="")); cat("\n")
}
sourceDir <- function(path, trace = FALSE, ...) {
for (nm in list.files(path, pattern = "\\.[RrSsQq]$")) {
if(trace) cat(nm,":")
source(file.path(path, nm), ...)
if(trace) cat("\n")
}
}
integrate <- function(f, a, b, n=1000) {
l = a+ ((b-a)/n)*0.5
A = 0
for (i in 0:(n-1)) {
A = A + f(l)*((b-a)/n)
l = l + ((b-a)/n)
}
return (A)
}
v <- function(i, n, moment=1) {
delta = log(1+i)*moment
return(exp(-delta*n))
}
tpx <- function(t,age=x) {
if( (age+t)>= w) return(0)
return(exp(-(A*t + B/log(c)*(c^age)*(c^t -1))))
}
tpselectx <- function(t,age=x) {
if( (age+t)>= w) return(0)
if (0<= t & t<=2) {
return(exp(0.9^(2-t)* ((A*(1 - 0.9^t)/log(0.9)) + (B*c^age*(c^t - 0.9^t)/log(0.9/c)))))
} else {
return(tpx(t,age))
}
}
tpselectxplus1 <- function(t,age=x) {
if( (age+t)>= w) return(0)
if(0<=t & t<=1) {
return(exp(0.9^(1-t)*( (A/log(0.9)*(1-0.9^t)) + (B*c^(age+1)/log(0.9/c)*(c^t - 0.9^t)))))
} else {
return(tpx(t,age))
}
}
tqx <- function(t,age=x) {
return(1-tpx(t,age))
}
tqselectx <- function(t,age=x) {
return(1-tpselectx(t,age))
}
tqselectxplus1 <- function(t,age=x) {
return(1-tpselectxplus1(t,age))
}
kdeferredqx <- function(k,age=x) {
return(tpx(k,age) - tpx(k+1,age))
}
kdeferredselectqx <- function(k,age=x) {
return(tpselectx(k,age) - tpselectx(k+1,age))
}
kdeferredselectqxplus1 <- function(k,age=x) {
return(tpselectxplus1(k,age)-tpselectxplus1(k+1,age))
}
uxplust <- function(t,age=x){
return(A+(B*c^(age+t)))
}
uxpluss <- function(t,age=x) {
if(t<0 || t>2) return(uxplust(t,age))
else return(0.9^(2-t)*uxplust(t,age))
}
udeferredtqx <- function(u,t,age=x) {
return(tpx(u,age) - tpx(u+t,age))
}
fxt <- function(t,age=x) {
return(tpx(t,age)*uxplust(t,age))
}
Axcolumn <- function(x,w,i,moment=1) {
axcol <- numeric(0); temp=0
for (k in 0:(w-x)) {
if(k==0) temp = v(i,1,moment)
else temp = axcol[k]*tpx(1,(w-k))*v(i,1,moment) + v(i,1,moment)*tqx(1,(w-k))
axcol = c(axcol, temp)
}
return(rev(axcol))
}
Aselectxcolumn <- function(x,w,i,moment=1) {
aselectxcol <- numeric(0); temp = 0
axplus2col = Axcolumn(x+2, w+2, i, moment)
for (k in 0:(w-x)) {
temp = (v(i,1,moment)*tqselectx(1,x+k) + v(i, 2, moment)*tqselectxplus1(1,x+k)*tpselectx(1,x+k)
+ v(i,2,moment)*tpselectx(2,x+k)*axplus2col[k+1])
aselectxcol = c(aselectxcol, temp)
}
return(aselectxcol)
}
aduexcolumn <- function(x,w,i,moment=1) {
aduexcol <- numeric(0); temp=0
axcol = Axcolumn(x,w,i,moment)
for (k in 0:(w-x)) {
aduexcol = c(aduexcol, (1-axcol[k])/d(i))
}
aduexcol = c(aduexcol,0)
return(aduexcol)
}
adueselectxcolumn <- function(x,w,i,moment=1) {
aduexcol <- numeric(0); temp=0
aselectxcol = Aselectxcolumn(x,w,i,moment)
for (k in 0:(w-x)) {
aduexcol = c(aduexcol, (1-aselectxcol[k])/d(i))
}
aduexcol = c(aduexcol,0)
return(aduexcol)
}
nEx <- function(n,x,moment=1) {
return(v(i,n,moment)*tpx(n,x))
}
nEselectx <- function(n,x,moment=1) {
return(v(i,n,moment)*tpselectx(n,x))
}
d <- function(i) {
return(i/(1+i))
}
iupperm <- function(i,m) {
return(m*((1+i)^(1/m) - 1))
}
TableD.1 <- function(x, w) {
age <- numeric(0)
lselectx <- numeric(0)
lselectxplus1 <- numeric(0)
lxplus2 <- numeric(0)
ageplus2 <- numeric(0)
for (j in (x-2) : (w-2) ) {
lxplus2 = c(lxplus2, l20*tpx( j-(x-2), x) )
ageplus2 = c(ageplus2, j+2)
if(j<x) {
age=c(age, 0)
lselectx=c(lselectx, 0)
lselectxplus1 = c(lselectxplus1, 0)
} else {
age = c(age,j)
lselectx = c(lselectx, lxplus2[j-(x-3)]/tpselectx(2,j) )
lselectxplus1 = c(lselectxplus1, lxplus2[j-(x-3)]/tpselectxplus1(1,j))
}
}
table <- data.frame(age, lselectx, lselectxplus1, lxplus2, ageplus2)
table = formatDataFrame(table, 2)
names(table) = c("x", "l[x]", "l[x]+1", "lx+2", "x+2")
return(table[0:(101-x),])
}
TableD.2 <- function(x,w) {
age <- numeric(0)
aselectxcol = Aselectxcolumn(x,w,i,1)
aselectx2col = Aselectxcolumn(x,w,i,2)
adueselectxcol = adueselectxcolumn(x,w,i,1)
Eselectx5 <- numeric(0)
Eselectx10 <- numeric(0)
Eselectx20 <- numeric(0)
for (j in (x:w)) {
age = c(age,j)
Eselectx5 = c(Eselectx5, nEselectx(5,j))
Eselectx10 = c(Eselectx10, nEselectx(10,j))
Eselectx20 = c(Eselectx20, nEselectx(20,j))
}
table <- data.frame(age, adueselectxcol, aselectxcol, aselectx2col, Eselectx5, Eselectx10, Eselectx20,age)
table = formatDataFrame(table, 5)
names(table) = c("x", '\"ax', "A[x]", "2A[x]", "5E[x]", "10E[x]", "20E[x]","x")
return(table[0:(81-x),])
}
TableD.3 <- function(x,w) {
age <- numeric(0)
aduexcol = aduexcolumn(x,w,i,1)
axcol = Axcolumn(x,w,i,1)
ax2col = Axcolumn(x,w,i,2)
Ex5 <- numeric(0)
Ex10 <- numeric(0)
Ex20 <- numeric(0)
for (j in (x:w)) {
age = c(age,j)
Ex5 = c(Ex5, nEx(5,j))
Ex10 = c(Ex10, nEx(10,j))
Ex20 = c(Ex20, nEx(20,j))
}
table <- data.frame(age, aduexcol, axcol, ax2col, Ex5, Ex10, Ex20,age)
table = formatDataFrame(table,5)
names(table) = c("x", '\"ax' ,"Ax", "2Ax", "5Ex", "10Ex", "20Ex","x")
return(table[0:(81-x),])
}
Axterm <- function(x,n, bar=FALSE, m=1, moment=1) {
# continuous case
if (bar == TRUE) {
# function to integrate
intFunction <- function(t) {
return(v(i,t,moment)*fxt(t,x))
}
return(integrate(intFunction,0,n,1000))
}
# mthly case
if (m>1) {
a <<- x
total = 0
for (k in (0:((m*n)-1))) {
total = total + v(i, ((k+1)/m),moment)*udeferredtqx(k/m,1/m,x)
}
return(total)
}
# annual case (m=1)
total = 0
for (k in (0:(n-1))) {
total = total + v(i,k+1,moment)*kdeferredqx(k,x)
}
return(total)
}
annuityDue <- function(i,n,bar=FALSE) {
# continuous case
if(bar==TRUE) {
# delta = ln(1+i)
return( (1-v(i,n))/log(1+i))
}
# mthly case
if(m>1) {
return( (1-v(i,n))/iupperm(i,m))
}
return( (1-v(i,n))/d(i))
}
annuityImmediate <- function(i,n) {
return(v(i,1)*annuityDue(i,n))
}
aduex <- function(x,bar=FALSE,m=1,moment=1) {
# continuous case
if(bar==TRUE) {
return((1-Ax(x,bar,m,moment))/log(1+i))
}
# annual case
return((1-Ax(x,bar,m,moment))/d(i))
}
ax <- function(x,bar=FALSE,m=1,moment=1) {
return(aduex(x,bar=FALSE,m=1,moment=1)-1)
}
Ax <- function(x,bar=FALSE, m=1, moment=1) {
return(Axterm(x,w,bar,m,moment))
}
ndeferredAx <- function(x, n, bar=FALSE, m=1, moment=1) {
return(Ax(x+n,bar,m,moment)*nEx(n,x))
}
Axendow <- function(x, n, bar=FALSE, m=1, moment=1, UDD=FALSE) {
if(UDD==TRUE) return(nEx(n,x,moment) + Axterm(x,n,bar,m,moment)*sanglen(i))
else return(nEx(n,x,moment)+ Axterm(x,n,bar,m,moment))
}
Aselectx <- function(x,moment=1) {
return(v(i,1,moment)*tqselectx(1,x) + v(i,2,moment)*tqselectxplus1(1,x)*tpselectx(1,x) +
v(i,2,moment)*tpselectx(2,x)*(Ax(x+2,moment)))
}
tableD.1 = TableD.1(x,w)
tableD.2 = TableD.2(x,w)
tableD.3 = TableD.3(x,w)
print(tableD.1)
x l[x] l[x]+1 lx+2 x+2
1 0 0 0 100000 20
2 0 0 0 99975 21
3 20 100000 99974 99950 22
4 21 99970 99948 99924 23
5 22 99945 99923 99898 24
6 23 99919 99896 99871 25
7 24 99893 99870 99844 26
8 25 99866 99842 99816 27
9 26 99838 99814 99787 28
10 27 99810 99786 99758 29
11 28 99781 99756 99727 30
12 29 99752 99726 99696 31
13 30 99721 99694 99663 32
14 31 99689 99661 99629 33
15 32 99656 99627 99594 34
16 33 99622 99592 99557 35
17 34 99586 99555 99518 36
18 35 99549 99516 99477 37
19 36 99510 99475 99433 38
20 37 99468 99431 99387 39
21 38 99424 99385 99338 40
22 39 99378 99336 99286 41
23 40 99328 99283 99230 42
24 41 99275 99227 99169 43
25 42 99218 99166 99104 44
26 43 99156 99101 99034 45
27 44 99090 99030 98958 46
28 45 99019 98953 98875 47
29 46 98941 98870 98784 48
30 47 98856 98779 98685 49
31 48 98764 98679 98576 50
32 49 98663 98570 98457 51
33 50 98553 98451 98326 52
34 51 98431 98319 98182 53
35 52 98297 98174 98022 54
36 53 98150 98014 97846 55
37 54 97987 97836 97651 56
38 55 97807 97640 97435 57
39 56 97608 97423 97196 58
40 57 97387 97182 96930 59
41 58 97142 96915 96634 60
42 59 96870 96618 96306 61
43 60 96568 96287 95941 62
44 61 96232 95920 95534 63
45 62 95859 95512 95083 64
46 63 95444 95057 94580 65
47 64 94981 94552 94020 66
48 65 94467 93989 93398 67
49 66 93895 93363 92706 68
50 67 93259 92668 91937 69
51 68 92551 91894 91082 70
52 69 91765 91035 90134 71
53 70 90891 90081 89082 72
54 71 89922 89024 87917 73
55 72 88847 87852 86628 74
56 73 87656 86556 85203 75
57 74 86340 85124 83633 76
58 75 84885 83546 81904 77
59 76 83283 81809 80006 78
60 77 81519 79901 77927 79
61 78 79584 77812 75657 80
62 79 77466 75532 73186 81
63 80 75154 73050 70507 82
64 81 72640 70360 67615 83
65 82 69916 67456 64507 84
66 83 66978 64337 61185 85
67 84 63825 61004 57657 86
68 85 60459 57465 53935 87
69 86 56889 53733 50039 88
70 87 53128 49829 45996 89
71 88 49198 45779 41841 90
72 89 45129 41620 37619 91
73 90 40956 37395 33380 92
74 91 36725 33157 29184 93
75 92 32491 28965 25094 94
76 93 28312 24883 21178 95
77 94 24253 20978 17502 96
78 95 20382 17316 14126 97
79 96 16764 13957 11103 98
80 97 13458 10954 8470 99
81 98 10515 8342 6248 100
print(tableD.2)
x "ax A[x] 2A[x] 5E[x] 10E[x] 20E[x] x
1 20 19.967 0.04918 0.00576 0.7825 0.6122 0.37440 20
2 21 19.921 0.05140 0.00610 0.7825 0.6122 0.37429 21
3 22 19.872 0.05373 0.00648 0.7825 0.6121 0.37417 22
4 23 19.820 0.05618 0.00689 0.7824 0.6121 0.37404 23
5 24 19.766 0.05874 0.00734 0.7824 0.6120 0.37390 24
6 25 19.710 0.06143 0.00783 0.7824 0.6120 0.37373 25
7 26 19.651 0.06424 0.00837 0.7824 0.6119 0.37354 26
8 27 19.589 0.06720 0.00895 0.7823 0.6118 0.37334 27
9 28 19.524 0.07029 0.00959 0.7823 0.6117 0.37310 28
10 29 19.456 0.07353 0.01028 0.7822 0.6116 0.37284 29
11 30 19.384 0.07693 0.01104 0.7822 0.6115 0.37254 30
12 31 19.310 0.08049 0.01186 0.7821 0.6114 0.37221 31
13 32 19.232 0.08421 0.01276 0.7821 0.6112 0.37183 32
14 33 19.150 0.08811 0.01373 0.7820 0.6111 0.37141 33
15 34 19.064 0.09220 0.01479 0.7819 0.6109 0.37094 34
16 35 18.974 0.09647 0.01594 0.7818 0.6107 0.37041 35
17 36 18.880 0.10094 0.01720 0.7817 0.6105 0.36982 36
18 37 18.782 0.10562 0.01856 0.7816 0.6102 0.36915 37
19 38 18.679 0.11051 0.02004 0.7814 0.6099 0.36841 38
20 39 18.572 0.11563 0.02164 0.7813 0.6096 0.36757 39
21 40 18.460 0.12097 0.02338 0.7811 0.6092 0.36663 40
22 41 18.342 0.12656 0.02527 0.7809 0.6088 0.36558 41
23 42 18.220 0.13240 0.02731 0.7807 0.6083 0.36440 42
24 43 18.092 0.13849 0.02952 0.7805 0.6078 0.36307 43
25 44 17.958 0.14485 0.03191 0.7802 0.6072 0.36159 44
26 45 17.819 0.15149 0.03450 0.7799 0.6066 0.35994 45
27 46 17.673 0.15841 0.03730 0.7796 0.6058 0.35809 46
28 47 17.522 0.16563 0.04032 0.7792 0.6050 0.35601 47
29 48 17.364 0.17314 0.04358 0.7788 0.6040 0.35370 48
30 49 17.200 0.18098 0.04709 0.7783 0.6030 0.35112 49
31 50 17.028 0.18913 0.05087 0.7777 0.6018 0.34824 50
32 51 16.850 0.19761 0.05495 0.7771 0.6005 0.34503 51
33 52 16.665 0.20642 0.05933 0.7764 0.5990 0.34146 52
34 53 16.473 0.21558 0.06404 0.7757 0.5974 0.33749 53
35 54 16.273 0.22509 0.06909 0.7748 0.5955 0.33308 54
36 55 16.066 0.23496 0.07451 0.7738 0.5934 0.32819 55
37 56 15.851 0.24519 0.08031 0.7727 0.5911 0.32279 56
38 57 15.628 0.25579 0.08653 0.7715 0.5885 0.31681 57
39 58 15.398 0.26677 0.09317 0.7701 0.5856 0.31024 58
40 59 15.160 0.27811 0.10025 0.7686 0.5823 0.30300 59
41 60 14.913 0.28984 0.10781 0.7669 0.5786 0.29508 60
42 61 14.659 0.30194 0.11586 0.7649 0.5746 0.28641 61
43 62 14.397 0.31442 0.12441 0.7628 0.5700 0.27698 62
44 63 14.127 0.32727 0.13350 0.7603 0.5650 0.26674 63
45 64 13.850 0.34049 0.14313 0.7576 0.5593 0.25569 64
46 65 13.564 0.35407 0.15333 0.7546 0.5531 0.24381 65
47 66 13.272 0.36801 0.16411 0.7511 0.5461 0.23112 66
48 67 12.972 0.38230 0.17548 0.7473 0.5384 0.21764 67
49 68 12.665 0.39692 0.18746 0.7430 0.5298 0.20343 68
50 69 12.351 0.41186 0.20005 0.7383 0.5204 0.18856 69
51 70 12.031 0.42710 0.21326 0.7329 0.5099 0.17313 70
52 71 11.705 0.44262 0.22709 0.7270 0.4985 0.15730 71
53 72 11.373 0.45840 0.24154 0.7204 0.4859 0.14122 72
54 73 11.037 0.47442 0.25662 0.7130 0.4722 0.12511 73
55 74 10.696 0.49065 0.27229 0.7048 0.4572 0.10918 74
56 75 10.352 0.50706 0.28856 0.6957 0.4409 0.09368 75
57 76 10.004 0.52362 0.30541 0.6857 0.4232 0.07887 76
58 77 9.654 0.54029 0.32280 0.6745 0.4043 0.06500 77
59 78 9.302 0.55704 0.34072 0.6622 0.3840 0.05230 78
60 79 8.950 0.57382 0.35912 0.6486 0.3624 0.04096 79
61 80 8.597 0.59061 0.37797 0.6337 0.3395 0.03113 80
print(tableD.3)
x "ax Ax 2Ax 5Ex 10Ex 20Ex x
1 20 19.966 0.04922 0.00580 0.7825 0.6122 0.37440 20
2 21 19.920 0.05144 0.00614 0.7825 0.6122 0.37429 21
3 22 19.871 0.05378 0.00652 0.7825 0.6121 0.37417 22
4 23 19.819 0.05622 0.00694 0.7824 0.6121 0.37404 23
5 24 19.765 0.05879 0.00739 0.7824 0.6120 0.37390 24
6 25 19.709 0.06147 0.00788 0.7824 0.6120 0.37373 25
7 26 19.650 0.06429 0.00841 0.7824 0.6119 0.37354 26
8 27 19.588 0.06725 0.00900 0.7823 0.6118 0.37334 27
9 28 19.523 0.07034 0.00964 0.7823 0.6117 0.37310 28
10 29 19.455 0.07359 0.01033 0.7822 0.6116 0.37284 29
11 30 19.383 0.07698 0.01109 0.7822 0.6115 0.37254 30
12 31 19.309 0.08054 0.01192 0.7821 0.6114 0.37221 31
13 32 19.230 0.08427 0.01281 0.7821 0.6112 0.37183 32
14 33 19.148 0.08817 0.01379 0.7820 0.6111 0.37141 33
15 34 19.063 0.09226 0.01486 0.7819 0.6109 0.37094 34
16 35 18.973 0.09653 0.01601 0.7818 0.6107 0.37041 35
17 36 18.879 0.10101 0.01727 0.7817 0.6105 0.36982 36
18 37 18.780 0.10569 0.01863 0.7816 0.6102 0.36915 37
19 38 18.678 0.11059 0.02012 0.7814 0.6099 0.36841 38
20 39 18.570 0.11571 0.02173 0.7813 0.6096 0.36757 39
21 40 18.458 0.12106 0.02347 0.7811 0.6092 0.36663 40
22 41 18.340 0.12665 0.02536 0.7809 0.6088 0.36558 41
23 42 18.218 0.13249 0.02741 0.7807 0.6083 0.36440 42
24 43 18.090 0.13859 0.02963 0.7805 0.6078 0.36307 43
25 44 17.956 0.14496 0.03203 0.7802 0.6072 0.36159 44
26 45 17.816 0.15161 0.03463 0.7799 0.6066 0.35994 45
27 46 17.671 0.15854 0.03744 0.7796 0.6058 0.35809 46
28 47 17.519 0.16577 0.04047 0.7792 0.6050 0.35601 47
29 48 17.361 0.17330 0.04374 0.7788 0.6040 0.35370 48
30 49 17.196 0.18114 0.04727 0.7783 0.6030 0.35112 49
31 50 17.025 0.18931 0.05108 0.7777 0.6018 0.34824 50
32 51 16.846 0.19780 0.05517 0.7771 0.6005 0.34503 51
33 52 16.661 0.20664 0.05957 0.7764 0.5990 0.34146 52
34 53 16.468 0.21582 0.06430 0.7757 0.5974 0.33749 53
35 54 16.268 0.22535 0.06938 0.7748 0.5955 0.33308 54
36 55 16.060 0.23524 0.07483 0.7738 0.5934 0.32819 55
37 56 15.844 0.24550 0.08067 0.7727 0.5911 0.32279 56
38 57 15.621 0.25613 0.08692 0.7715 0.5885 0.31681 57
39 58 15.390 0.26714 0.09360 0.7701 0.5856 0.31024 58
40 59 15.151 0.27852 0.10073 0.7686 0.5823 0.30300 59
41 60 14.904 0.29028 0.10834 0.7669 0.5786 0.29508 60
42 61 14.649 0.30243 0.11644 0.7649 0.5746 0.28641 61
43 62 14.386 0.31495 0.12506 0.7628 0.5700 0.27698 62
44 63 14.115 0.32785 0.13421 0.7603 0.5650 0.26674 63
45 64 13.836 0.34113 0.14392 0.7576 0.5593 0.25569 64
46 65 13.550 0.35477 0.15420 0.7546 0.5531 0.24381 65
47 66 13.256 0.36878 0.16507 0.7511 0.5461 0.23112 66
48 67 12.954 0.38313 0.17654 0.7473 0.5384 0.21764 67
49 68 12.646 0.39783 0.18862 0.7430 0.5298 0.20343 68
50 69 12.330 0.41285 0.20133 0.7383 0.5204 0.18856 69
51 70 12.008 0.42818 0.21467 0.7329 0.5099 0.17313 70
52 71 11.680 0.44379 0.22864 0.7270 0.4985 0.15730 71
53 72 11.347 0.45968 0.24324 0.7204 0.4859 0.14122 72
54 73 11.008 0.47580 0.25847 0.7130 0.4722 0.12511 73
55 74 10.665 0.49215 0.27433 0.7048 0.4572 0.10918 74
56 75 10.318 0.50868 0.29079 0.6957 0.4409 0.09368 75
57 76 9.967 0.52536 0.30783 0.6857 0.4232 0.07887 76
58 77 9.614 0.54217 0.32544 0.6745 0.4043 0.06500 77
59 78 9.260 0.55906 0.34359 0.6622 0.3840 0.05230 78
60 79 8.904 0.57599 0.36224 0.6486 0.3624 0.04096 79
61 80 8.548 0.59293 0.38134 0.6337 0.3395 0.03113 80
#### TESTING OF PROGRAM ####
# ages = c(20,40,60,80);col1<-numeric(0);col2<-numeric(0);col3<-numeric(0)
# for(age in ages) {
# col1 = c(col1, Axterm(age,10,bar=TRUE))
# col2 = c(col2, Axterm(age,10,m=4))
# col3 = c(col3, Axterm(age,10))
# }
# example4.4 = data.frame(col1,col2,col3)
# print(example4.4)
#
# #### exercise 4.2 ####
# a = Axterm(30,20)
# b = Axendow(40,20,UDD=TRUE)
# c = ndeferredAx(25,10)