-
Notifications
You must be signed in to change notification settings - Fork 2
/
utility_functions.R
186 lines (153 loc) · 6.36 KB
/
utility_functions.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
library(zoo)
#Present Value function (for constant payments)
pv = function(rate, nper, pmt) {
PV = (1 + rate)*pmt * (1 - (1 + rate) ^ (-nper)) / rate
return(PV)
}
#NPV function
npv = function(rate, cashflows) {
for(i in 1:length(cashflows)){
if(i == 1){
NPV <- cashflows[i]/((1+rate)^(i))
} else {
NPV <- NPV + cashflows[i]/((1+rate)^(i))
}
}
return(NPV)
}
#Amo payment functions
#pmt0 = basic amo payment calculation, assuming payment beginning of period
PMT0 <- function(r, nper, pv) {
if (r == 0) {
a <- pv/nper
} else {
a <- pv*r*(1+r)^(nper-1)/((1+r)^nper-1)
}
# if(nper == 0){
# a <- 0
# }
return(a)
}
#pmt = amo payment function with growth rate and timing added; t = 1 for end of period payment, 0.5 for half period, and 0 for beginning of period.
PMT <- function(r, g = 0, nper, pv, t = 1) {
a <- PMT0((1+r)/(1+g) - 1, nper, pv*(1+r)^t)
return(a)
}
#Funding period function
NPER <- function(r,g,pv,t,pmt){
PV <- pv*(1+r)^t
R <- (1+r)/(1+g) - 1
TempValue <- PV*R/(pmt*(1+R))
NPER <- -log(1-TempValue,(1+R))
if(is.infinite(NPER)) {NPER <- 100}
if(is.nan(NPER)) {NPER <- 100}
return(NPER)
}
# NPER(0.07, g = 0, pv = 5000, t = 1, pmt = 1000)
#Cumulative future values function
cumFV <- function(interest, cashflow){
cumvalue <- double(length = length(cashflow))
for (i in 2:length(cumvalue)) {
cumvalue[i] <- cumvalue[i - 1]*(1 + interest) + cashflow[i - 1]
}
return(cumvalue)
}
#Geometric average return function
geo_return <- function(x, na.rm = F) {
if (na.rm == T) {
x = na.omit(x)
}
avg_return <- prod(1+x)^(1/length(x)) - 1
return(avg_return)
}
#Estimate geometric average return from arithmetic return and standard deviation (annualized)
est_geo_return <- function(mean_return, sd_return){
geo_return <- mean_return - sd_return^2/(2*(1+mean_return))
return(geo_return)
}
#Estimate arithmetic return from geometric return and standard deviation (annualized)
root_geo_return <- function(mean_return, sd_return, geometric_return){
est_geo_return(mean_return, sd_return) - geometric_return
}
est_arith_return <- function(geometric_return, sd_return){
x <- uniroot(f = root_geo_return,
interval = c(-1,1),
geometric_return = geometric_return,
sd_return = sd_return)
return(x$root)
}
#Function to determine the interest crediting rate of a cash balance plan
smooth_return <- function(returns, floor, cap, upside_share) {
x <- min(max(floor, floor + upside_share * (geo_return(returns) - floor)), cap)
return(x)
}
expected_icr <- function(geometric_return, sd_return, smooth_period, floor, cap, upside_share, n_periods = 30, n_simulations = 10000) {
mean_return <- est_arith_return(geometric_return, sd_return)
simulated_returns <- matrix(rnorm(n_simulations * n_periods, mean = mean_return, sd = sd_return), nrow = n_periods, ncol = n_simulations)
initial_returns <- matrix(floor, nrow = smooth_period - 1, ncol = n_simulations)
returns_matrix <- rbind(initial_returns, simulated_returns)
smooth_returns <- rollapply(returns_matrix, width = smooth_period, align = "right", FUN = smooth_return, floor = floor, cap = cap, upside_share = upside_share)
avg_interest_rates <- apply(smooth_returns, 2, FUN = geo_return)
exp_interest_rate <- median(avg_interest_rates)
return(exp_interest_rate)
}
# expected_icr(geometric_return = 0.08, sd_return = 0.16, smooth_period = 5, floor = 0.06, cap = 0.1, upside_share = 0.7, n_periods = 30)
#Adding new entrants function
add_new_entrants <- function(g, ne_dist, wf1, wf2, ea, age, position_matrix){
#g is the assumed population growth of the plan
#ne_dist is a vector representing the distribution of new entrants for each entry age
#wf1 is the population in period 1. wf2 is the wf1 population after decremented
#ea and age are two vectors representing entry age and age for active members
#position_matrix is the matrix that rearranges the new entrant numbers in the right positions to be added to the active workforce array
ne <- sum(wf1)*(1 + g) - sum(wf2)
ne_vec <- ne * ne_dist
ne_matrix <- matrix(ne_vec, nrow = length(ea), ncol = length(age))
ne_matrix_trans <- ne_matrix * position_matrix
return(ne_matrix_trans)
}
#Recursive growing function (with lag)
recur_grow <- function(x, g) {
if (length(x) > 1) {
for (i in 2:length(x)) {
x[i] <- x[i-1] * (1 + g[i - 1])
}
}
return(x)
}
#Recursive growing function (no lag)
recur_grow2 <- function(x, g) {
if (length(x) > 1) {
for (i in 2:length(x)) {
x[i] <- x[i-1] * (1 + g[i])
}
}
return(x)
}
#Present Value of Future Benefits (PVFB) function (to be applied to a vector of "Pension Wealth") for active members
#sep_rate_vec is a vector containing separation rates. Interest is a single discount rate (ARR). value_vect is a vector containing the present values of pension benefits at separation ages.
#The purpose of this function is to calculate the PVFB at each active age (not just the entry age)
PVFB <- function(sep_rate_vec, interest, value_vec) {
PVFB <- double(length = length(value_vec))
for (i in 1:length(value_vec)) {
sep_rate <- sep_rate_vec[i:length(sep_rate_vec)]
#sep_prob in a given year is the probability that the member will survive all the previous years and get terminated exactly in the given year
sep_prob <- cumprod(1 - lag(sep_rate, n = 2, default = 0)) * lag(sep_rate, default = 0)
value <- value_vec[i:length(value_vec)]
value_adjusted <- value * sep_prob
PVFB[i] <- npv(interest, value_adjusted[2:length(value_adjusted)])
}
return(PVFB)
}
#Present Value of Future Salaries (PVFS) function (to be applied to a vector of salaries)
#remaining_prob_vec is a vector containing the remaining probabilities. Interest is a si ngle discount rate (ARR). sal_vec is a vector containing the salaries.
PVFS <- function(remaining_prob_vec, interest, sal_vec) {
PVFS <- double(length = length(sal_vec))
for (i in 1:length(sal_vec)) {
remaining_prob_og <- remaining_prob_vec[i:length(remaining_prob_vec)]
remaining_prob <- remaining_prob_og / remaining_prob_og[1]
sal <- sal_vec[i:length(sal_vec)]
sal_adjusted <- sal * remaining_prob
PVFS[i] <- npv(interest, sal_adjusted)
}
return(PVFS)
}