Retirement : simulating wealth with random returns, inflation and withdrawals – Shiny web application

Today, I want to share the Retirement : simulating wealth with random returns, inflation and withdrawals – Shiny web application (code at GitHub).

This application was developed and contributed by Pierre Chretien, I only made minor updates. This is application is a great example of how easy it is to convert your R script into interactive Shiny web application.

Please see below the sample script to simulate wealth random returns, inflation and withdrawals:

#-------------------------------------
# Inputs
#-------------------------------------

# Initial capital
start.capital = 2000000

# Investment
annual.mean.return = 5 / 100
annual.ret.std.dev = 7 / 100

# Inflation
annual.inflation = 2.5 / 100
annual.inf.std.dev = 1.5 / 100

# Withdrawals
monthly.withdrawals = 10000

# Number of observations (in Years)
n.obs = 20

# Number of simulations
n.sim = 200

#-------------------------------------
# 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 Withdrawals
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.withdrawals
}

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

# convert to millions
nav = nav / 1000000

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

# plot all scenarios    
matplot(nav, type = 'l', las = 1, xlab = 'Months', ylab = 'Millions', 
	main = 'Projected Value of initial capital')

# plot % of scenarios that are still paying
p.alive = 1 - rowSums(is.na(nav)) / n.sim

plot(100 * p.alive, las = 1, xlab = 'Months', ylab = 'Percentage Paying', 
	main = 'Percentage of Paying Scenarios', ylim=c(0,100))
grid()	

# plot distribution of final wealth
final.nav = nav[n.obs + 1, ]
	final.nav = final.nav[!is.na(final.nav)]

plot(density(final.nav, from=0, to=max(final.nav)), las = 1, xlab = 'Final Capital', 
	main = paste('Distribution of Final Capital,', 100 * p.alive[n.obs + 1], '% are still paying'))
grid()	

The corresponding Shiny web application consists of two files:

• ui.r – User Interface
• server.r – R simulations and calculations

Following is the user interface (ui.r) that maps and describes required inputs for the retirement simulation:

# Define UI for application that plots random distributions 
shinyUI(pageWithSidebar(

  # Application title
  headerPanel("Retirement : simulating wealth with random returns, inflation and withdrawals"),

  # Sidebar with a slider input for number of observations
  sidebarPanel(
    sliderInput("n.obs", 
                "Number of observations (in Years):", 
                min = 0, 
                max = 40, 
                value = 20),

    sliderInput("start.capital", 
                "Initial capital invested :", 
                min = 100000, 
                max = 10000000, 
                value = 2000000,
                step = 100000,
                format="$#,##0",
                locale="us"),

    sliderInput("annual.mean.return", 
                "Annual return from investments (in %):", 
                min = 0.0, 
                max = 30.0, 
                value = 5.0,
                step = 0.5),

    sliderInput("annual.ret.std.dev", 
                "Annual volatility from investments (in %):", 
                min = 0.0, 
                max = 25.0, 
                value = 7.0, 
                step = 0.1),

    sliderInput("annual.inflation", 
                "Annual inflation (in %):", 
                min = 0, 
                max = 20, 
                value = 2.5,
                step = 0.1),

    sliderInput("annual.inf.std.dev", 
                "Annual inflation volatility. (in %):", 
                min = 0.0, 
                max = 5.0,
                value = 1.5,
                step = 0.05),

    sliderInput("monthly.withdrawals", 
                "Monthly capital withdrawals:", 
                min = 1000, 
                max = 100000, 
                value = 10000,
                step = 1000,
                format="$#,##0",
                locale="us",),
                
    sliderInput("n.sim", 
                "Number of simulations:", 
                min = 0, 
                max = 2000, 
                value = 200)
                
  ),

  # Show a plot of the generated distribution
  mainPanel(
    plotOutput("distPlot", height = "600px")
  )
))

The last step is modify the retirement simulation logic to use user inputs:

library(shiny)

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

  # Function that generates scenarios and computes NAV.
  getNav <- reactive({ 
	#-------------------------------------
	# Inputs
	#-------------------------------------
	
	# Initial capital
	start.capital = input$start.capital
	
	# Investment
	annual.mean.return = input$annual.mean.return / 100
	annual.ret.std.dev = input$annual.ret.std.dev / 100
	
	# Inflation
	annual.inflation = input$annual.inflation / 100
	annual.inf.std.dev = input$annual.inf.std.dev / 100
	
	# Withdrawals
	monthly.withdrawals = input$monthly.withdrawals
	
	# Number of observations (in Years)
	n.obs = input$n.obs
	
	# Number of simulations
	n.sim = input$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 Withdrawals
	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.withdrawals
	}	
	
	# once nav is below 0 => run out of money
	nav[ nav < 0 ] = NA
	
	# convert to millions
	nav = nav / 1000000
	
	return(nav)  
  })
  
  # Expression that plot NAV paths. 
  output$distPlot <- renderPlot({
	nav = getNav()

	layout(matrix(c(1,2,1,3),2,2))
	
	# plot all scenarios    
	matplot(nav, type = 'l', las = 1, xlab = 'Months', ylab = 'Millions', 
		main = 'Projected Value of initial capital')

		
	# plot % of scenarios that are still paying
	p.alive = 1 - rowSums(is.na(nav)) / ncol(nav)
	
	plot(100 * p.alive, las = 1, xlab = 'Months', ylab = 'Percentage Paying', 
		main = 'Percentage of Paying Scenarios', ylim=c(0,100))
	grid()	

		
	last.period = nrow(nav)
  	
	# plot distribution of final wealth
	final.nav = nav[last.period, ]
		final.nav = final.nav[!is.na(final.nav)]
	
	if(length(final.nav) ==  0) return()		
	
	plot(density(final.nav, from=0, to=max(final.nav)), las = 1, xlab = 'Final Capital', 
		main = paste('Distribution of Final Capital,', 100 * p.alive[last.period], '% are still paying'))
	grid()	
  })
		
})

We all done now!!! Shiny is amazing in the way it allows you to convert your script into interactive web application with just two simple steps.

Please play around with the Retirement : simulating wealth with random returns, inflation and withdrawals – Shiny web application (code at GitHub).

Have a good weekend

Maximum Sharpe Portfolio

Maximum Sharpe Portfolio or Tangency Portfolio is a portfolio on the efficient frontier at the point where line drawn from the point (0, risk-free rate) is tangent to the efficient frontier.

There is a great discussion about Maximum Sharpe Portfolio or Tangency Portfolio at quadprog optimization question. In general case, finding the Maximum Sharpe Portfolio requires a non-linear solver because we want to find portfolio weights w to maximize w' mu / sqrt( w' V w ) (i.e. Sharpe Ratio is a non-linear function of w). But as discussed in the quadprog optimization question, there are special cases when we can use quadratic solver to find Maximum Sharpe Portfolio. As long as all constraints are homogeneous of degree 0 (i.e. if we multiply w by a number, the constraint is unchanged – for example, w1 > 0 is equivalent to 5*w1 > 5*0), the quadratic solver can be used to find Maximum Sharpe Portfolio or Tangency Portfolio.

I have implemented the logic to find Maximum Sharpe Portfolio or Tangency Portfolio in the max.sharpe.portfolio() function in strategy.r at github. You can specify following 2 parameters:

• Type of portfolio: ‘long-only’, ‘long-short’, or ‘market-neutral’
• Portfolio exposure. For example, ‘long-only’ with exposure = 1, is a fully invested portfolio

Now, let’s construct a sample efficient frontier and plot Maximum Sharpe Portfolio.

###############################################################################
# Load Systematic Investor Toolbox (SIT)
# http://systematicinvestor.wordpress.com/systematic-investor-toolbox/
###############################################################################
setInternet2(TRUE)
con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb'))
    source(con)
close(con)
 
	#*****************************************************************
	# Create Efficient Frontier
	#****************************************************************** 	
	# create sample historical input assumptions
	ia = aa.test.create.ia()
	
	# create long-only, fully invested efficient frontier
	n = ia$n		

	# 0 <= x.i <= 1
	constraints = new.constraints(n, lb = 0, ub = 1)
		constraints = add.constraints(diag(n), type='>=', b=0, constraints)
		constraints = add.constraints(diag(n), type='<=', b=1, constraints)
		
	# SUM x.i = 1
	constraints = add.constraints(rep(1, n), 1, type = '=', constraints)		
	
	# create efficient frontier
	ef = portopt(ia, constraints, 50, 'Efficient Frontier') 
	
	#*****************************************************************
	# Create Plot
	#****************************************************************** 	
	# plot efficient frontier
	plot.ef(ia, list(ef), transition.map=F)	 
	
	# find maximum sharpe portfolio
	max(portfolio.return(ef$weight,ia) /  portfolio.risk(ef$weight,ia))
	
	# plot minimum variance portfolio
	weight = min.var.portfolio(ia,constraints)	
	points(100 * portfolio.risk(weight,ia), 100 * portfolio.return(weight,ia), pch=15, col='red')
	portfolio.return(weight,ia) /  portfolio.risk(weight,ia)
		
	# plot maximum Sharpe or tangency portfolio
	weight = max.sharpe.portfolio()(ia,constraints)	
	points(100 * portfolio.risk(weight,ia), 100 * portfolio.return(weight,ia), pch=15, col='orange')
	portfolio.return(weight,ia) /  portfolio.risk(weight,ia)
		
	plota.legend('Minimum Variance,Maximum Sharpe','red,orange', x='topright')	

Now let’s see how to construct ‘long-only’, ‘long-short’, or ‘market-neutral’ Maximum Sharpe Portfolio or Tangency Portfolios:

	#*****************************************************************
	# Examples of Maximum Sharpe or Tangency portfolios construction
	#****************************************************************** 	
	weight = max.sharpe.portfolio('long-only')(ia,constraints)	
		round(weight,2)
		round(c(sum(weight[weight<0]), sum(weight[weight>0])),2)
		
	weight = max.sharpe.portfolio('long-short')(ia,constraints)			
		round(weight,2)
		round(c(sum(weight[weight<0]), sum(weight[weight>0])),2)
		
	weight = max.sharpe.portfolio('market-neutral')(ia,constraints)			
		round(weight,2)
		round(c(sum(weight[weight<0]), sum(weight[weight>0])),2)	

The long-only Maximum Sharpe portfolio as expected has exposure of 100%. The long-short Maximum Sharpe portfolio is 227% long and 127% short. The market-neutral Maximum Sharpe portfolio is 100% long and 100% short.

As the last step, I run Maximum Sharpe algo vs other portfolio optimization methods I have previously discussed (i.e. Risk Parity, Minimum Variance, Maximum Diversification, Minimum Correlation) on the 10 asset universe used in the Adaptive Asset Allocation post.

	#*****************************************************************
	# Load historical data
	#****************************************************************** 
	load.packages('quantmod')
	
	tickers = spl('SPY,EFA,EWJ,EEM,IYR,RWX,IEF,TLT,DBC,GLD')

	data <- new.env()
	getSymbols(tickers, src = 'yahoo', from = '1980-01-01', env = data, auto.assign = T)
		for(i in ls(data)) data[[i]] = adjustOHLC(data[[i]], use.Adjusted=T)							
	bt.prep(data, align='keep.all', dates='2004:12::')
 
	#*****************************************************************
	# Code Strategies
	#******************************************************************
	prices = data$prices  
	n = ncol(prices)
   
	models = list()
   
 	#*****************************************************************
	# Code Strategies
	#******************************************************************
	# find period ends
	period.ends = endpoints(prices, 'months')
        period.ends = period.ends[period.ends > 0]
        
	n.mom = 180
	n.vol = 60
	n.top = 4        
	momentum = prices / mlag(prices, n.mom)  
       
	obj = portfolio.allocation.helper(data$prices, period.ends=period.ends,
		lookback.len = n.vol, universe = ntop(momentum[period.ends,], n.top) > 0,
		min.risk.fns = list(EW=equal.weight.portfolio,
						RP=risk.parity.portfolio,
						MV=min.var.portfolio,
						MD=max.div.portfolio,
						MC=min.corr.portfolio,
						MC2=min.corr2.portfolio,
						MCE=min.corr.excel.portfolio,
						MS=max.sharpe.portfolio())
	) 
	
	models = create.strategies(obj, data)$models
					
	#*****************************************************************
	# Create Report
	#******************************************************************    
	strategy.performance.snapshoot(models, T)

	plotbt.custom.report.part2(models$MS)

	# Plot Portfolio Turnover for each strategy
	layout(1)
	barplot.with.labels(sapply(models, compute.turnover, data), 'Average Annual Portfolio Turnover')

The allocation using Maximum Sharpe Portfolio produced more concentrated portfolios with higher total return, higher Sharpe ratio, and higher turnover.

Cluster Risk Parity back-test

In the Cluster Portfolio Allocation post, I have outlined the 3 steps to construct Cluster Risk Parity portfolio. At each rebalancing period:

• Create Clusters
• Allocate funds within each Cluster using Risk Parity
• Allocate funds across all Clusters using Risk Parity

I created a helper function distribute.weights() function in strategy.r at github to automate these steps. It has 2 parameters:

• Function to allocate funds. For example, risk.parity.portfolio, will use use risk parity to allocate funds both within and across clusters.
• Function to create clusters. For example, cluster.group.kmeans.90, will create clusters using k-means algorithm

Here is the example how to put it all together. Let’s first load historical prices for the 10 major asset classes:

###############################################################################
# Load Systematic Investor Toolbox (SIT)
# http://systematicinvestor.wordpress.com/systematic-investor-toolbox/
###############################################################################
setInternet2(TRUE)
con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb'))
    source(con)
close(con)

	#*****************************************************************
	# Load historical data for ETFs
	#****************************************************************** 
	load.packages('quantmod')

	tickers = spl('GLD,UUP,SPY,QQQ,IWM,EEM,EFA,IYR,USO,TLT')

	data <- new.env()
	getSymbols(tickers, src = 'yahoo', from = '1900-01-01', env = data, auto.assign = T)
		for(i in ls(data)) data[[i]] = adjustOHLC(data[[i]], use.Adjusted=T)
		
	bt.prep(data, align='remove.na')

Next, let’s run the 2 versions of Cluster Portfolio Allocation using Equal Weight and Risk Parity algorithms to allocate funds:

	#*****************************************************************
	# Code Strategies
	#****************************************************************** 	
	periodicity = 'months'
	lookback.len = 250
	cluster.group = cluster.group.kmeans.90
	
	obj = portfolio.allocation.helper(data$prices, 
		periodicity = periodicity, lookback.len = lookback.len,
		min.risk.fns = list(
				EW=equal.weight.portfolio,
				RP=risk.parity.portfolio,
						
				C.EW = distribute.weights(equal.weight.portfolio, cluster.group),
				C.RP=distribute.weights(risk.parity.portfolio, cluster.group)
			)
	) 		
	
	models = create.strategies(obj, data)$models

Finally, let’s examine the results:

	#*****************************************************************
	# Create Report
	#****************************************************************** 	
	strategy.performance.snapshoot(models, T)

The Cluster Portfolio Allocation produce portfolios with better risk-adjusted returns and smaller drawdowns.

To view the complete source code for this example, please have a look at the bt.cluster.portfolio.allocation.test1() function in bt.test.r at github.

Sector Rotation Back Test Shiny web application

Today, I want to share the Sector Rotation Back Test application (code at GitHub).

This is the last application in the series of examples (I have shared 5 examples) that will demonstrate the amazing Shiny framework and Systematic Investor Toolbox to analyze stocks, make back-tests, and create summary reports.

The motivation for this series of posts is to show how to translate, with minimal effort, your back-test scripts into live web applications that can either be run from the Shiny server or your personal computer.

Please note that Back Test applications take longer time to update the plots/tables and hence maybe more appropriate to run from your local computer.

For example to run Sector Rotation application from your computer, you first need to install Shiny. And next execute following commands:

library(shiny)
runGitHub('SIT','systematicinvestor', subdir = 'Shiny/sector.rotation')

Market Filter Back Test Shiny web application

Today, I want to share the Market Filter Back Test application (code at GitHub).

This is the forth application in the series of examples (I plan to share 5 examples) that will demonstrate the amazing Shiny framework and Systematic Investor Toolbox to analyze stocks, make back-tests, and create summary reports.

The motivation for this series of posts is to show how to translate, with minimal effort, your back-test scripts into live web applications that can either be run from the Shiny server or your personal computer.

Please note that Back Test applications take longer time to update the plots/tables and hence maybe more appropriate to run from your local computer.

I will post the next example application on Monday.

Have a good long weekend