R Optimization Returning Incorrect Values

Question

My objective function:

helper.post<-function(monthly.mean.return,
                  start.capital,  #initial nest egg
                  target.legacy,
                  monthly.inflation.post,
                  monthly.withdrawals,
                  n.obs){ 

  req = matrix(start.capital, n.obs+1, 1) #matrix for storing target weight

  for (a in 1:n.obs) {
    #cat("a: ",a,"\n")
    req[a + 1, ] = req[a, ] * (1 + monthly.mean.return - monthly.inflation.post) -     monthly.withdrawals[a,]
  }
  ending.value=req[nrow(req),]
  #ending.value
  value=target.legacy - ending.value

  return(abs(value))
}

With the following Optimization structure, changing the n.obs between the two values give the same output:

ie if n.obs = 288 or n.obs = 336, it gives the same values.

optimize(f=helper.post,
     start.capital = 1000000,
     target.legacy = 1000000,
     monthly.inflation.post=0.002083333,
     monthly.withdrawals = matrix(rep(10000,n.obs)),
     n.obs = n.obs, 
     lower = 0,
     upper = 1,
     tol = 0.00000000000000000000000000000000001)$minimum

The value is correct seems to be a estimation as oppose to the correct value. Any idea what I may be doing incorrectly? Would a different optimization tool work better for such precise optimization efforts? I tried uni-root, but it doesn't sem to work as the end points are not opposite signs..

uniroot( helper.post, 
     c(0, 1),
     start.capital = start.capital,
     target.legacy = target.legacy,
     monthly.inflation.post=monthly.inflation.post,
     monthly.withdrawals = monthly.withdrawals,
     n.obs = n.obs)$root

Answer 1

Let's start with a slight rewrite of your code. I replaced one-column matrices with vectors. I also added an option for returning the error itself or its absolute value. You'll want to use the absolute value when trying to minimize the error with optim while you'll want the value itself when trying to find its root with uniroot.

helper.post <- function(monthly.mean.return,
                        start.capital,
                        target.legacy,
                        monthly.inflation.post,
                        monthly.withdrawals,
                        n.obs,
                        return.abs = TRUE) { 

  req <- numeric(n.obs + 1)
  req[1] <- start.capital

  for (month in 1:n.obs) {
    req[month + 1] <- req[month] *
                      (1 + monthly.mean.return - monthly.inflation.post) -
                      monthly.withdrawals[month]
  }
  ending.value <- req[n.obs + 1]
  error <- target.legacy - ending.value

  return(ifelse(return.abs, abs(error), error))
}

Now let's optimize it:

n.obs <- 288

optimize(f = helper.post,
         start.capital = 1000000,
         target.legacy = 1000000,
         monthly.inflation.post = 0.002083333,
         monthly.withdrawals = matrix(rep(10000,n.obs)),
         n.obs = n.obs, 
         lower = 0,
         upper = 1,
         tol = 1e-20)$minimum
# [1] 0.01208333

And let's check the solution with uni.root:

uniroot(helper.post, 
        c(0, 1),
        start.capital = 1000000,
        target.legacy = 1000000,
        monthly.inflation.post = 0.002083333,
        monthly.withdrawals = matrix(rep(10000,n.obs)),
        n.obs = n.obs,
        return.abs = FALSE,
        tol = 1e-20)$root
# [1] 0.01208333

They match. There is nothing wrong with one or the other tool...

If you run again with a different value for n.obs, you will get the exact same result. Why? Because you have picked constant withdrawals and equal values for the start and target capitals: the output you get is the rate needed to maintain that balance constant from one month to the next, regardless of the total number of months.

In fact, this is the kind of thing you would want to put in a unit test. Because it is an expected and easily interpretable result. Another one that comes to mind is if you made the withdrawals equal to zero everywhere. Then you would expect your answer to match the inflation rate. Give it a try and see that it does indeed.