In grfiv/sp500SlidingWindow: Sliding Window Investment Analysis

knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(collapse = TRUE,              # hadley
                      comment  = "#>",              # hadley
                      error=TRUE, purl=FALSE,       # to be able to see errors
                      fig.width=7.25, fig.height=6) # nice-sized pictures

Scenario

We are facing retirement, having to live off a lump sum we have invested. How much can we withdraw each year adjusted for inflation?

Parameters

* we are 65 years old   
* we have a nest egg of $3,000,000
* we expect to live for 30 years   
* we expect 2.5% annual inflation            (average since 1914 = 3.27%)
* we expect to pay brokerage fees of 1.25%   (US average per Morningstar) 
* Treasuries have a 2.16% coupon             (as of May 2016)
* our all-in ordinary-income tax rate is 20%

# function-specific parameters
library(sp500SlidingWindow)
window_width <- 30           # life expectancy
annual_fee   <- 0.0125       # brokerage fees
output_path  <- paste0(tempdir(), "/")

# analysis-specific parameters
current_age        <- 65            # starting age
annual_inflation   <- 0.025         # annual inflation
initial_investment <- 3000000       # our nest egg

# withdraw this much each year
withdrawal_percent <- 0.03171
initial_withdrawal <- initial_investment * withdrawal_percent   

# for bond analysis
coupon   <- 0.0216
tax_rate <- 0.20

Investment vector

We retire with this much money to live on for the rest of our 'window_width'-year life

investment_vector <- c(initial_investment,
                       rep(0,window_width-1))

Withdrawal vector

We hope to spend the initial_withdrawal adjusted for inflation each year

withdrawal_vector <- sapply(0:(window_width-1), function(i) {
    return(initial_withdrawal * (1 + annual_inflation)**i)
})

# total amount we hope to withdraw over the whole period
total_hoped_for_wdr <- sum(withdrawal_vector)

# print year-by-year hoped-for withdrawals
knitr::kable(data.frame(age        = seq(current_age, current_age+window_width-1),  
                        withdrawal = fmt(withdrawal_vector)))
# save as image
library(gridExtra)
png(filename=paste0(output_path, "wdr_vector.png"), width=250, height=700)
  grid.table(data.frame(age        = seq(current_age, current_age+window_width-1),  
                        withdrawal = fmt(withdrawal_vector)))
dev.off()

Total withdrawals over r paste0(window_width, " years: $", fmt(total_hoped_for_wdr))

Keep the Money in Cash

If we do not invest the lump sum but simply spend it down, what will be the result?

remaining_amt    <- initial_investment
ages             <- current_age:(window_width+current_age-1)
remaining_amount <- sapply(1:window_width, function(i) {
    remaining_amt <<- remaining_amt - withdrawal_vector[i]
    return(remaining_amt)
})

# test if we go below zero
if (length(which(remaining_amount<=0))) {
    # goes below zero
    plot_points = TRUE
    busted_age <- current_age+which(remaining_amount<=0)[1]-1
    sub = paste0("Run out of money at age ", busted_age)
} else {
    # does not go below zero
    plot_points = FALSE
    sub="Starting amount covers withdrawals without investment"
}

#print(remaining_amount)

ylim <- c(min(min(remaining_amount, 0)), initial_investment)
plot(ages, remaining_amount, pch=20,
     xlab="Age", ylab=NA, yaxt="n", ylim=ylim,
     main="What if I simply live off the cash?",
     sub=sub)
if (plot_points) {
    points(which(remaining_amount<=0)+current_age-1, 
           remaining_amount[which(remaining_amount<=0)], col="red", pch=20)
    abline(v=busted_age, col="red", lty=2)
    abline(h=0, col="red", lty=2)
}

axis(2, las=2, at=axTicks(2), labels=fmt(axTicks(2)), cex.axis=0.72)
grid()

dev.copy(png, paste0(output_path, "invested_in_cash.png"))
dev.off()

Buy Bonds

Why not just buy a US Treasury bond, pay the taxes on the interest payments and draw down on the remainder? (Somewhat simplified because you can't sell part of a single bond, you have to buy several bonds so you can liquidate part of the portfolio when you need to.)

bond_df <- buy_bonds(investment_vector,
                     withdrawal_vector,
                     coupon,
                     tax_rate)

#knitr::kable(bond_df)

pos_wd <- bond_df$bal[which(bond_df$bal > 0)]

# test if we go below zero
if (length(pos_wd) < nrow(bond_df)) {
    # goes below zero
    plot_points   <- TRUE
    busted_age    <- current_age+which(bond_df$bal <= 0)[1]-1
    sub           <- paste0("Run out of money at age ", busted_age)
    remaining_amt <- cumsum(-withdrawal_vector[(length(pos_wd)+1):length(withdrawal_vector)])
    ylim          <- c(min(remaining_amt), initial_investment)
} else {
    # does not go below zero
    plot_points = FALSE
    sub="Starting amount covers withdrawals without investment"
    ylim          <- c(0, initial_investment)
}

plot(ages, bond_df$bal, type="n", pch=20,
     xlab="Age", ylab=NA, yaxt="n", ylim=ylim,
     main="What if I buy a bond?",
     sub=sub)

axis(2, las=2, at=axTicks(2), labels=fmt(axTicks(2)), cex.axis=0.72)
grid()

points(ages[1:length(pos_wd)], pos_wd, pch=20)

if (plot_points) {
    points(ages[(length(pos_wd)+1):nrow(bond_df)], remaining_amt, pch=20, col="red")
    abline(h=0,          lty=2, col="red")
    abline(v=busted_age, lty=2, col="red")
}

dev.copy(png, paste0(output_path, "invest_in_a_bond.png"))
dev.off()

Invest in the Stock Market

The idea of sliding-window analysis is to ask how a certain set of annual investments and withdrawals would perform in each of the periods of a certain width of the stock market.

We are simulating a brokerage cash account so the balance can never go below zero. If the withdrawal vector calls for an amount greater than the current balance, the remaining balance is withdrawn and the balance is set to zero.

The analysis in this example invests the lump sum at the beginning of each window and tracks the effect of the stock market on that investment, making inflation-adjusted withdrawals each year. The critical question is whether the investor will run out of money in any of the periods we test.

window_df <- sp500SlidingWindow(investment_vector,
                                withdrawal_vector,
                                window_width = window_width,
                                annual_fee   = annual_fee,
                                output_path  = output_path)

#knitr::kable(window_df)

Number of failed periods r length(which(window_df$ending_bal <= 0)); r paste0(round((1-(length(which(window_df$ending_bal <= 0))/nrow(window_df)))*100,0), "%") percent success

How are the results distributed?

plot_ending_bal_distribution(window_df, NULL, window_width)

Worst Period

library(png)
worst_period <- which.min(window_df$wdr)

worst_year_path <- paste0(output_path, 
                          window_df$start_year[worst_period], "-",
                          window_df$end_year[worst_period], ".png")
pp <- readPNG(worst_year_path, native = TRUE, info = TRUE)
plot(0:1, 0:1, type="n", ann=FALSE, axes=FALSE)
rasterImage(pp, 0, 0, 1, 1)

All Failed Periods

library(dplyr)
failed_df     <- window_df %>% filter(ending_bal <= 0) %>% arrange(wdr)
failed_df$inv <- fmt(failed_df$inv)
failed_df$wdr <- fmt(failed_df$wdr)
failed_df$IRR <- ifelse(is.na(failed_df$IRR), NA, paste0(round(failed_df$IRR*100, 2), "%"))
knitr::kable(failed_df)

Best Period

best_period  <- which.max(window_df$ending_bal)

best_year_path <- paste0(output_path, 
                          window_df$start_year[best_period], "-",
                          window_df$end_year[best_period], ".png")
pp <- readPNG(best_year_path, native = TRUE, info = TRUE)
plot(0:1, 0:1, type="n", ann=FALSE, axes=FALSE)
rasterImage(pp, 0, 0, 1, 1)

Modal (Most-Likely) Period

bal_hist_obj    <- hist(window_df$ending_bal, plot=FALSE)
mod_end_bal     <- bal_hist_obj$mids[which.max(bal_hist_obj$counts)]
mod_window_idx  <- which.min(abs(window_df$ending_bal - mod_end_bal))
closest_mod_bal <- window_df$ending_bal[mod_window_idx]

mode_year_path <- paste0(output_path, 
                          window_df$start_year[mod_window_idx], "-",
                          window_df$end_year[mod_window_idx], ".png")
pp <- readPNG(mode_year_path, native = TRUE, info = TRUE)
plot(0:1, 0:1, type="n", ann=FALSE, axes=FALSE)
rasterImage(pp, 0, 0, 1, 1)

Middling Periods

middling_df     <- window_df %>%  arrange(ending_bal, wdr, start_year)
lmid            <- which.min(abs(middling_df$ending_bal - closest_mod_bal))
middling_df     <- middling_df[(lmid-2):(lmid+2),]
middling_df$inv <- fmt(middling_df$inv)
middling_df$wdr <- fmt(middling_df$wdr)
middling_df$IRR <- ifelse(is.na(middling_df$IRR), NA, paste0(round(middling_df$IRR*100, 2), "%"))
middling_df$ending_bal <- fmt(middling_df$ending_bal)
knitr::kable(middling_df)

Inflation

library(rvest)

inflation_table <- read_html("http://www.usinflationcalculator.com/inflation/historical-inflation-rates/") %>%
    html_node("table") %>%
  html_table(header = TRUE)

for (col in names(inflation_table)) {
    inflation_table[, col] <- as.numeric(inflation_table[, col])
}

Average annual inflation rate since r inflation_table$Year[1]: r paste0(round(mean(inflation_table$Ave, na.rm = TRUE), 2),"%")

Environment

devtools::session_info()

grfiv/sp500SlidingWindow documentation built on May 17, 2019, 8:36 a.m.