401K simulator

401K Simulator

You go into work day after day and you have some money taken out every week/month whatever. How is it growing? Is it going to be enough to retire when you’re 90 years old? Did I make a mistake about saying 90 years, didn’t I really mean 62? Well, in these uncertain times we live in (when were the certain times?) you’ll want to know how that hard earned cash is growing, or at least how much it should be growing. That’s where this tool comes into play. You put in a bunch of settings and it spits out a bunch of simulations of what could happen. Since your money is probably put into some dodgy 401K fund that your company so lovingly provides, you’re gonna want to know what is the pay off in X years from now. Or in other words, when can i realistically ditch this job that is sucking the life out of my soul. Play around with some of the variables and find out how long it will take you to cash in and head to the beach for good. Or, unless you’re one of those people that loves their job and does not want to retire, you might want to just get back to work since you love it so much. What are you doing here!

Calculations and initial parameters for simulating 401K growth.

  • Years to maturity, retirement
  • Initial Capital
  •  Annual rate of return
  •  Annual risk
  •  Inflation volatily
  •  Monthly deposits
  •  Number of simulations

This simulation is an adaptation of the retirement app from Systematic Investor.
Here is the R source code for the server.R file:

library(shiny)

paramNames <- c("start_capital", "annual_mean_return", "annual_ret_std_dev",
"annual_inflation", "annual_inf_std_dev", "monthly_deposits", "n_obs",
"n_sim")

# Define server logic required to generate and plot a random distribution
#
# Idea and original code by Pierre Chretien
# Small updates by Michael Kapler
# modified by eddy o'donnell
#
shinyServer(function(input, output, session) {

getParams <- function(prefix) {
input[[paste0(prefix, "_recalc")]]

params <- lapply(paramNames, function(p) {
input[[paste0(prefix, "_", p)]]
})
names(params) <- paramNames
params
}

  # Function that generates scenarios and computes NAV. The expression
  # is wrapped in a call to reactive to indicate that:
  #
  # 1) It is "reactive" and therefore should be automatically
  # re-executed when inputs change
  #
  navA <- reactive(do.call(simulate_nav, getParams("a")))

  # Expression that plot NAV paths. The expression
  # is wrapped in a call to renderPlot to indicate that:
  #
  # 1) It is "reactive" and therefore should be automatically
  # re-executed when inputs change
  # 2) Its output type is a plot
  #
  output$a_distPlot <- renderPlot({
   plot_nav(navA(), output)
  })
  

})

simulate_nav <- function(start_capital = 2000000, annual_mean_return = 5.0,
                         annual_ret_std_dev = 7.0, annual_inflation = 2.5,
                         annual_inf_std_dev = 1.5, monthly_deposits = 1000,
                         n_obs = 20, n_sim = 200) {
  #-------------------------------------
  # Inputs convert the percentables by dividing by 100
  #-------------------------------------

  # get percentage - Investment
  annual.mean.return = annual_mean_return / 100
  annual.ret.std.dev = annual_ret_std_dev / 100

  # get percentage -  Inflation
  annual.inflation = annual_inflation / 100
  annual.inf.std.dev = annual_inf_std_dev / 100

  # deposits = monthly_deposits
  # Number of observations (in Years) = n.obs
  # Number of simulations = n_sim

  #-------------------------------------
  # Simulation
  #-------------------------------------

  # number of months to simulate
  n_obs = 12 * n_obs

  # monthly Investment and Inflation assumptions
  monthly.mean.return = annual.mean.return / 12
  monthly.ret.std.dev = annual.ret.std.dev / sqrt(12)

  monthly.inflation = annual.inflation / 12
  monthly.inf.std.dev = annual.inf.std.dev / sqrt(12)

  # simulate Returns
  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)

  # simulate deposits
  nav = matrix(start_capital, n_obs + 1, n_sim)
  for (j in 1:n_obs) {
    nav[j + 1, ] = nav[j, ] * (1 + monthly.invest.returns[j, ] - monthly.inflation.returns[j, ]) + monthly_deposits
  }

  # once nav is below 0 => run out of money
  nav[ nav < 0 ] = NA

  # convert to thousands
  nav = nav / 1000

  return(nav)
}

plot_nav <- function(nav, output) {

  # layout(matrix(c(1,2,1,3),2,2))

  # palette(c("red", "grey50", "grey30", "grey70", "#d9230f"))
  now <- as.numeric(format(Sys.time(),"%Y"))
  n <- nrow(nav) 
  yearLabels <- seq(now+1,now + round(n/12),1)
  l <- length(yearLabels)
  top <-round(n/l)
  myat <- seq(10,n,top)

  # plot all scenarios
  matplot(nav,
    type = 'l', lwd = 0.5, lty = 1, col = rainbow(ncol(nav)),
    xaxt = 'n', ylab = 'Dollar Amount in Thousands', sub='Projected Value Over Time', cex.sub=2)
    axis(3,at=myat,labels=yearLabels,cex.axis=1.25,tck=0.01)
    grid()


  last.period = nrow(nav)

  # plot distribution of final wealth
  final.nav = nav[last.period, ]
  final.nav = final.nav[!is.na(final.nav)]
  final.nav = final.nav * 1000
  if(length(final.nav) == 0) return()
  #h<-hist(final.nav, col='#2470CF',xlab="Dollar amount in thousands.",main="Histogram of the last period of the simulation")
  # xfit<-seq(min(final.nav),max(final.nav),length=length(final.nav))
  # yfit<-dnorm(xfit,mean=mean(final.nav),sd=sd(final.nav))
  # yfit <- yfit*diff(h$mids[1:2])*length(final.nav)
  #lines(xfit, yfit, col='#202020', lwd=2)

  #plot(density(final.nav, from=0, to=max(final.nav)), las = 1, xlab = 'Final Capital',
  #  main = paste0('Distribution of Final Capital\n', 100 * p.alive[last.period], '% are still paying'))
  #grid()
  output$summary <- renderPrint({
    summary(final.nav)
  })

}

Leave a Reply

Your email address will not be published. Required fields are marked *

2 × 5 =