CODEWITHSUNDEEP
roptimization
user1234440
user1234440

Reputation: 23607

R Optimization given objective function

obj1<-function(monthly.savings,
               success,
               start.capital,
               target.savings,
               monthly.mean.return,
               monthly.ret.std.dev,
               monthly.inflation,
               monthly.inf.std.dev,
               n.obs,
               n.sim=1000){

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

  monthly.invest.returns = matrix(0, n.obs, n.sim)
  monthly.inflation.returns = matrix(0, n.obs, n.sim)

  monthly.invest.returns[] = rnorm(n.obs * n.sim, mean = monthly.mean.return, sd = monthly.ret.std.dev)
  monthly.inflation.returns[] = rnorm(n.obs * n.sim, mean = monthly.inflation, sd = monthly.inf.std.dev)

  #for loop to be 
  for (a in 1:n.obs){
    req[a + 1, ] = req[a, ] * (1 + monthly.invest.returns[a,] - monthly.inflation.returns[a,]) + monthly.savings
  }

  ending.values=req[nrow(req),]
  suc<-sum(ending.values>target.savings)/n.sim
  value<-success-suc

  return(abs(value))
}

I have the above objective function that I want to minimize for. It tries to solve for the monthly savings required for a given probability of success. Given the following input assumptions

success<-0.9
start.capital<-1000000
target.savings<-1749665
monthly.savings=10000
monthly.mean.return<-(5/100)/12
monthly.ret.std.dev<-(3/100)/sqrt(12)
monthly.inflation<-(5/100)/12
monthly.inf.std.dev<-(1.5/100)/sqrt(12)
monthly.withdrawals<-10000
n.obs<-10*12  #years * 12 months in a year
n.sim=1000

I used the following notation:

optimize(f=obj1,
        success=success,
        start.capital=start.capital,
        target.savings=target.savings,
        monthly.mean.return=monthly.mean.return,
        monthly.ret.std.dev=monthly.ret.std.dev,
        monthly.inflation=monthly.inflation,
        monthly.inf.std.dev=monthly.inf.std.dev,
        n.obs = n.obs,
        n.sim = n.sim,
        lower = 0,
        upper = 10000,
        tol = 0.000000001,maximum=F)

I get 7875.03

Since I am sampling from a normal distribution, the output will be different each time but they should be around the same give or take a few % points. The problem I am having is that I can't specify a upper limit arbitrarily. The above example's upper limit (10000) is cherry picked after numerous trials. If say I put in a upper limit of 100000 (unreasonable I know) it will return that number as oppose to finding the global minimum saving. Any ideas where I am structuring my objective function incorrectly?

thanks,

Upvotes: 0

Views: 1671

Answers (1)

Vincent Zoonekynd
Vincent Zoonekynd

Reputation: 32401

The fact that your function does not always return the same output for a given input is likely to pose a few problems (it will create a lot of spurious local minima): you can avoid them by setting the seed of the random number generator inside the function (e.g., set.seed(1)), or by storing the random numbers and reusing them each time, or by using a low-discrepancy sequence (e.g., randtoolbox::sobol).

Since it is a function of one variable, you can simply plot it to see what happens: it has a plateau after 10,000 -- optimization algorithms cannot distinguish between a plateau and a local optimum.

f <- function(x) {
  set.seed(1)
  obj1(x,
      success             = success,
      start.capital       = start.capital,
      target.savings      = target.savings,
      monthly.mean.return = monthly.mean.return,
      monthly.ret.std.dev = monthly.ret.std.dev,
      monthly.inflation   = monthly.inflation,
      monthly.inf.std.dev = monthly.inf.std.dev,
      n.obs               = n.obs,
      n.sim               = n.sim 
  )
}
g <- Vectorize(f)
curve(g(x), xlim=c(0, 20000))

Your initial problem is actually not a minimization problem, but a root finding problem, which is much easier.

obj2 <- function(monthly.savings) {
  set.seed(1)
  req = matrix(start.capital, n.obs+1, n.sim)
  monthly.invest.returns <- matrix(0, n.obs, n.sim)
  monthly.inflation.returns <- matrix(0, n.obs, n.sim)
  monthly.invest.returns[] <- rnorm(n.obs * n.sim, mean = monthly.mean.return, sd = monthly.ret.std.dev)
  monthly.inflation.returns[] <- rnorm(n.obs * n.sim, mean = monthly.inflation, sd = monthly.inf.std.dev)
  for (a in 1:n.obs)
    req[a + 1, ] <- req[a, ] * (1 + monthly.invest.returns[a,] - monthly.inflation.returns[a,]) + monthly.savings
  ending.values <- req[nrow(req),]
  suc <- sum(ending.values>target.savings)/n.sim
  success - suc
}
uniroot( obj2, c(0, 1e6) )
# [1] 7891.187

Upvotes: 4

Related Questions