I am trying to find a function input value (scalar) that, when combined with three other known vector parameters, creates a fourth vector that diminishes to zero exactly at a specific vector element (and not before).
The fourth vector contains calculated portfolio balances and each vector element is a year. I want the portfolio balance to diminish to zero exactly in the portfolio balance for the year the investor dies (the year is specified in one of the known vector parameters). I want to achieve this by finding the initial portfolio value that results in that vector element equalling zero in the year of death and no sooner.
(Note: In Excel, this could be done with the Goal Seek tool by setting a specific portfolio balance element (year) equal to zero by changing the initial portfolio value.)
I have tried both optim and optimize but last tried the script below. The optimized initial portfolio value always equals the constraint max.
Any help would be greatly appreciated. Been pulling out my hair for two days.
f1 <- function(init,spendV,returnV,stateV) {
# create a vector that equals an initial portfolio value (init) and sets subsequent vector elements
# equal to the previous portfolio balance minus an amount of spending and then applies a market return rate.
# Portfolio balance can never fall below zero
# When the state vector element = 4, the investor died in the previous period
endState <- which (stateV == 4)
portfoliobal <- rep(0,length(spendV))
portfoliobal[1] <- (init - spendV[1]) * returnV[1]
for (i in (2:length(spendV))) {
portfoliobal[i] <- max(0,(portfoliobal[i-1] - spendV[i]) * returnV[i])
}
# the function returns the portfolio balance when the investor dies
return(portfoliobal[endState])
}
# spending, market growth rates, and the period of death are given vectors
spendV <- c(97719.19,97737.92,102649.4,98669.15,78108.44,58105.49,51710.02,53267.12,
39982.34,21070,22439.68,21375.4,15613.45,10826.54,9333.82,7808.239,
8737.435,8382.001,7267.976,6534.688,5129.403,5026.947,4931.132,5468.401,
5033.245,5195.273,5199.938,4854.684,5039.221,3757.753,1822.97,1202.24,
1237.238,965.111,1051.235,906.4884,1110.66,1018.127,500.788,538.2703,
584.5545,599.1832,575.2254,640.8828,604.0179,781.878,595.9795,625.5037,
615.471,667.4227)
init <- 1000000
returnV <- c(0.9347388,1.170053,1.204515,0.9572276,1.044682,0.9229759,1.110595,0.9299398,
1.161509,1.053207,1.058104,0.8997761,1.000342,1.353597,1.031785,1.121795,0.8745584,
1.05637,1.180234,0.9795393,0.9137375,0.772738,1.021843,0.9697467,1.055284,0.9182615,
0.9662726,1.105152,1.099005,0.9195565,0.895424,0.9226368,1.196467,1.085768,0.9529325,
1.485245,0.9124764,0.8978044,0.8021779,0.9064698,1.034353,0.9914232,0.742632,
0.9308539,0.9683604,0.9325817,1.000051,1.145982,1.018012,1.127159)
stateV <- c(3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,
2,2,2,2,2,2,2,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
cat("\n\nCritical vector element= ",which(stateV == 4))
# the critical vector element is the portfolio balance at death
# I experiment with different inital portfolio values and find that portfolio balance[34] approaches zero
# when initial portfolio balance is about 710,000
f1value <- f1(init=1000000,spendV,returnV,stateV)
cat ("\n\nf1value for critical vector element",which(stateV == 4)," when inital value = 1,000,000: ",f1value)
f1value <- f1(init=800000,spendV,returnV,stateV)
cat ("\n\nf1value for critical vector element",which(stateV == 4)," when inital value = 800,000: ",f1value)
f1value <- f1(init=710000,spendV,returnV,stateV)
cat ("\n\nf1value for critical vector element",which(stateV == 4)," when inital value = 710,000: ",f1value)
# I am trying to find the initial portfolio balance that leaves zero dollars in portfolio balance vector [34]
# ---- BUT not BEFORE portfolio balance [34]
# opt <- optimize(f1,c(0,10000000),spendV,returnV,stateV,maximum = FALSE)
Answer
optim() and optimize() both seek a minimum (by default). You want a zero. The function you are looking for is uniroot(). Do this:
uniroot(f=f1, interval=c(1, 1e7), spendV=spendV, returnV=returnV, stateV=stateV)
#> $root
#> 703932.9
Some other stuff
For fun, here's a re-write of your function in vectorized notation (I lazily left out the returnV and spendV args):
f2 <- function(init, end=34) {
init*prod(returnV[1:end]) - sum(cumprod(returnV[end:1])*spendV[end:1])
}
uniroot(f2, c(1, 1e7))
#> $root
#> 703932.9
In general, if you are testing possible solutions, you might find it easier to plot a whole bunch of solutions, rather than trying them one at a time and reporting the result with cat(). Here's how I might go about it:
x <- seq(1, 2e6, by=100000)
y <- f2(x, 34)
plot(x, y, pch=20)
abline(v=0, h=0)
One last thought: in your function f1, the line which(stateV==4) might return a vector of values if the user inputs, say, stateV=c(1,4,4). I'm guessing you want the first occurrence of a "4", which you can accomplish with which(stateV==4)[1] or if(sum(stateV==4)>1) stop("Too many 4s!").