Archive

Archive for September, 2012

Weekend Reading – Gold in October

September 29, 2012 Leave a comment

I recently came across the “An early Halloween for gold traders” article by Mark Hulbert. I have discussed this type of seasonality analysis in my presentation at R/Finance this year.

It is very easy to run the seasonality analysis using the Systematic Investor Toolbox.

###############################################################################
# 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
    #****************************************************************** 
    load.packages('quantmod')
    ticker = 'GLD'
    
    data = getSymbols(ticker, src = 'yahoo', from = '1970-01-01', auto.assign = F)
        data = adjustOHLC(data, use.Adjusted=T)
        
    #*****************************************************************
    # Look at the Month of the Year Seasonality
    #****************************************************************** 
    month.year.seasonality(data, ticker)

This confirms that October have been historically bad for Gold, but we used only 8 years of history because GLD only started traded in 2004.

To get a more complete picture, there is a long history of Gold prices at the Bundes Bank. I found this data source used at the Wikiposit.

I created a helper function bundes.bank.data.gold() function in data.r at github to download prices from the Bundes Bank site.

    #*****************************************************************
    # Load long series of gold prices from Bundes Bank
    #****************************************************************** 
    data = bundes.bank.data.gold()

    #*****************************************************************
    # Look at the Month of the Year Seasonality
    #****************************************************************** 
    month.year.seasonality(data, 'GOLD', lookback.len = nrow(data))

The October have been historically bad for Gold using longer time series as well.

Next I would recommend looking at the daily Gold’s performance in October to get a better picture. You might want to use the Seasonality Tool for this purpose. Please read the Historical Seasonality Analysis: What company in DOW 30 is likely to do well in January? post for a case study on how to use the Seasonality Tool.

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

Categories: R

Calling Minimum Correlation Algorithm from Excel using RExcel & VBA

September 27, 2012 1 comment

I want to show the example of calling the Minimum Correlation Algorithm from Excel. I will use RExcel to connect R and Excel and will create a small VBA cell array function to communicate between Excel and R.

I have previously discussed the concept of connecting R and Excel in the “Calling Systematic Investor Toolbox from Excel using RExcel & VBA” post. Please read this post for the instructions to setup RExcel.

Following is a screen shot of the complete workbook:

You can download the MinimumCorrelation.xls workbook and experiment with it while you keep reading.

I created the “MinimumCorrelation” cell array function in VBA. Do not forget to use CTRL+SHIFT+ENTER to enter “MinimumCorrelation” function into your workbooks. The “MinimumCorrelation” function will send historical price information from Excel to the R environment. It will next execute the R script to construct weights using the Minimum Correlation Algorithm, and finally it will collect R calculations of portfolio weights and transfer them back to Excel.

Here is the R script that calls min.corr.portfolio() function. I created a VBA function “create_rcode” to create this file automatically for this example.

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

   #*****************************************************************          
   # Setup                                                                     
   #*****************************************************************          
   n = ncol(hist_prices)                                                       
   hist = na.omit(hist_prices / mlag(hist_prices) - 1)                         
                                                                               
   # create historical input assumptions                                       
   ia = list()                                                                 
       ia$n = n                                                                
       ia$risk = apply(hist, 2, sd)                                            
       ia$correlation = cor(hist, use='complete.obs', method='pearson')        
       ia$cov = ia$correlation * (ia$risk %*% t(ia$risk))                      

   # portfolio allocation                                                      
   weights = min.corr.portfolio(ia, null)                                      
   dim(weights)=c(1,n)                                                         

Next, the “MinimumCorrelation” cell array function in VBA that calls min.corr.portfolio() function in R:

'Minimum Correlation Algorithm
Public Function MinimumCorrelation(ByRef r_values As Range) As Variant
    ' Start R connection
    RInterface.StartRServer

    ' Write R code to file
    create_rcode

    ' Put Historical Asset Prices into R
    RInterface.PutArray "hist_prices", r_values
    
    ' Executes the commands in filename
    RInterface.RunRFile r_filename
         
    ' Get Portfolio Allocation determined by the Minimum Correlation Algorithm into Excel
    MinimumCorrelation = RInterface.GetArrayToVBA("weights")
End Function

The complete working copy of the MinimumCorrelation.xls workbook.

Please do not forget to use CTRL+SHIFT+ENTER to enter “MinimumCorrelation” function into your workbooks.

Categories: R

Minimum Correlation Algorithm Speed comparison

September 26, 2012 Leave a comment

The Minimum Correlation Algorithm is a heuristic method discovered by David Varadi. Below I will benchmark the execution speed of 2 versions of the Minimum Correlation Algorithm versus the traditional minimum variance optimization that relies on solving a quadratic programming problem.

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

	#*****************************************************************
	# Setup test input assumptions
	#*****************************************************************
	load.packages('quadprog,corpcor')
	
	n = 100
	hist = matrix(rnorm(1000*n), nc=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 historical input assumptions
	ia = list()
		ia$n = n
		ia$risk = apply(hist, 2, sd)
		ia$correlation = cor(hist, use='complete.obs', method='pearson')
		ia$cov = ia$correlation * (ia$risk %*% t(ia$risk))
				
		ia$cov = make.positive.definite(ia$cov, 0.000000001)
		ia$correlation = make.positive.definite(ia$correlation, 0.000000001)
		
	#*****************************************************************
	# Time each Algorithm 
	#*****************************************************************				
	load.packages('rbenchmark')			

	benchmark(
		min.var.portfolio(ia, constraints),
		min.corr.portfolio(ia, constraints),
		min.corr2.portfolio(ia, constraints),
		
		
	columns=c("test", "replications", "elapsed", "relative"),
	order="relative",
	replications=100
	)

I have run the code above for n=10 (10 assets), n=100 (100 assets), n=500 (500 assets), n=1000 (1000 assets)
[Please note that for n=1000 I have only run 5 replication]

n=10 (10 assets)
                      replications elapsed relative
   min.var.portfolio          100    0.02      1.0
 min.corr2.portfolio          100    0.02      1.0
  min.corr.portfolio          100    0.03      1.5


n=100 (100 assets)
                      replications elapsed relative
 min.corr2.portfolio          100    0.07 1.00
  min.corr.portfolio          100    0.25 3.57
   min.var.portfolio          100    0.31 4.42


n=500 (500 assets)
                      replications elapsed  relative
 min.corr2.portfolio          100    2.18  1.00
  min.corr.portfolio          100    9.59  4.39
   min.var.portfolio          100  139.13 63.82


n=1000 (1000 assets) 
                      replications elapsed relative
 min.corr2.portfolio            5    0.25     1.00
  min.corr.portfolio            5    1.39     5.56
   min.var.portfolio            5  113.27   453.08

For small universe (i.e. n ~ 100) all algorithms are fast. But once we attempt to solve 500 or 1000 assets portfolio allocation problem, the minimum variance algorithm is many times slower than the both versions of the minimum correlation algorithm.

So if we are considering a 500 assets weekly back-test for the 10 yrs the run-times in seconds (i.e. 52*10*single tun-time):

                            elapsed
 min.corr2.portfolio          11.3
  min.corr.portfolio          49.8
   min.var.portfolio         723.4

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

Categories: Portfolio Construction, R

Minimum Correlation Algorithm Example

September 24, 2012 4 comments

Today I want to follow up with the Minimum Correlation Algorithm Paper post and show how to incorporate the Minimum Correlation Algorithm into your portfolio construction work flow and also explain why I like the Minimum Correlation Algorithm.

First, let’s load the ETF’s data set used in the Minimum Correlation Algorithm Paper using the Systematic Investor Toolbox.

###############################################################################
# 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,quadprog')
	tickers = spl('SPY,QQQ,EEM,IWM,EFA,TLT,IYR,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='2002:08::')

Next, I created the portfolio.allocation.helper() function to ease the back-testing of portfolio construction algorithms:

	#*****************************************************************
	# Code Strategies
	#****************************************************************** 	
	
	obj = portfolio.allocation.helper(data$prices, periodicity = 'weeks',
		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),
		custom.stats.fn = 'portfolio.allocation.custom.stats'
	) 
	
	models = create.strategies(obj, data)$models

Please note that I assigned acronyms to various portfolio allocation algorithms in the code above.

For example, the Minimum Correlation Algorithm’s acronym is MC and the actual function that does all computations is min.corr.portfolio() function.

Next let’s create various summary reports to see the performance and allocations of different algorithms:

    #*****************************************************************
    # Create Report
    #******************************************************************       
	# Plot perfromance
	layout(1:2)
	plotbt(models, plotX = T, log = 'y', LeftMargin = 3)	    	
		mtext('Cumulative Performance', side = 2, line = 1)
		
	out = plotbt.strategy.sidebyside(models, return.table=T)

	# Plot time series of components of Composite Diversification Indicator
	cdi = custom.composite.diversification.indicator(obj,plot.table = F)	
		out = rbind(colMeans(cdi, na.rm=T), out)
		rownames(out)[1] = 'Composite Diversification Indicator(CDI)'
				
	# Portfolio Turnover for each strategy
	y = 100 * sapply(models, compute.turnover, data)
		out = rbind(y, out)
		rownames(out)[1] = 'Portfolio Turnover'		

	# Plot and compare strategies across different metrics
	performance.barchart.helper(out, 'Sharpe,Cagr,DVR,MaxDD', c(T,T,T,T))
	
	performance.barchart.helper(out, 'Volatility,Portfolio Turnover,Composite Diversification Indicator(CDI)', c(F,F,T))

	
	# Plot transition maps
	layout(1:len(models))
	for(m in names(models)) {
		plotbt.transition.map(models[[m]]$weight, name=m)
			legend('topright', legend = m, bty = 'n')
	}
	
	# Plot transition maps for Risk Contributions
	dates = index(data$prices)[obj$period.ends]
	for(m in names(models)) {
		plotbt.transition.map(make.xts(obj$risk.contributions[[m]], dates), 
		name=paste('Risk Contributions',m))
			legend('topright', legend = m, bty = 'n')
	}

Looking at the summary statistics, you mind wonder what makes the the Minimum Correlation Algorithm different and why do I like it?

The overall characteristics of the Minimum Correlation Algorithm makes it attractive for me. Specifically it has reasonable perfromance, limited draw-down, low turnover and high composite diversification score. The combination of these attractive properties does differentiate the Minimum Correlation Algorithm from others.

In the next post, I will show the speed benchmarks for the Minimum Correlation Algorithm.

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

Minimum Correlation Algorithm Paper

September 22, 2012 4 comments

Over summer I was busy collaborating with David Varadi on the Minimum Correlation Algorithm paper. Today I want to share the results of our collaboration:

The Minimum Correlation Algorithm is fast, robust, and easy to implement. Please add it to you portfolio construction toolbox and share your experience.

Categories: Portfolio Construction, R

Permanent Portfolio

September 18, 2012 2 comments

First, just a quick update: I’m moving the release date of the SIT package a few months down the road, probably in November.

Now back to the post. Recently I came across a series of interesting posts about the Permanent Portfolio at the GestaltU blog. Today I want to show you how to back-test the Permanent Portfolio using the Systematic Investor Toolbox.

The simple version of the Permanent Portfolio consists of equal allocations to stocks(SPY), gold(GLD), treasuries(TLT), and cash(SHY). [25% allocation each] The portfolio is rebalanced once a year if any allocation breaks out from the 15% – 35% range.

###############################################################################
# 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
	#****************************************************************** 
	load.packages('quantmod')	
	tickers = spl('SPY,TLT,GLD,SHY')
		
	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)
		
		# extend GLD with Gold.PM - London Gold afternoon fixing prices
		data$GLD = extend.GLD(data$GLD)
	
	bt.prep(data, align='remove.na')

	#*****************************************************************
	# Code Strategies
	#****************************************************************** 		
	prices = data$prices   
		n = ncol(prices)
		nperiods = nrow(prices)

	# annual
	period.ends = endpoints(prices, 'years')
		period.ends = period.ends[period.ends > 0]		
		period.ends.y = c(1, period.ends)

	# quarterly
	period.ends = endpoints(prices, 'quarters')
		period.ends = period.ends[period.ends > 0]		
		period.ends.q = c(1, period.ends)
					
	models = list()
	
	#*****************************************************************
	# Code Strategies
	#****************************************************************** 	
	target.allocation = matrix(rep(1/n,n), nrow=1)
	
	# Buy & Hold	
	data$weight[] = NA	
		data$weight[period.ends.y[1],] = target.allocation
	models$buy.hold = bt.run.share(data, clean.signal=F)
		
	# Equal Weight Annual
	data$weight[] = NA
		data$weight[period.ends.y,] = ntop(prices[period.ends.y,], n)
	models$equal.weight.y = bt.run.share(data, clean.signal=F)

	# Equal Weight Quarterly
	data$weight[] = NA
		data$weight[period.ends.q,] = ntop(prices[period.ends.q,], n)
	models$equal.weight.q = bt.run.share(data, clean.signal=F)

To rebalance base on the 10% threshold (i.e. portfolio weights breaking out from the 15% – 35% range) I will use bt.max.deviation.rebalancing() function introduced in the Backtesting Rebalancing methods post.

	#*****************************************************************
	# Rebalance based on threshold
	#****************************************************************** 	
	# Rebalance only when threshold is broken
	models$threshold.y = bt.max.deviation.rebalancing(data, models$buy.hold, target.allocation, 10/100, 0, period.ends = period.ends.y) 

	# Rebalance only when threshold is broken
	models$threshold.q = bt.max.deviation.rebalancing(data, models$buy.hold, target.allocation, 10/100, 0, period.ends = period.ends.q) 

	#*****************************************************************
	# Create Report
	#******************************************************************       
	plotbt.custom.report.part1(models)       
	
	plotbt.strategy.sidebyside(models)
		
	# Plot Portfolio Turnover for each Rebalancing method
	layout(1:2)
	barplot.with.labels(sapply(models, compute.turnover, data), 'Average Annual Portfolio Turnover', F)
	barplot.with.labels(sapply(models, compute.max.deviation, target.allocation), 'Maximum Deviation from Target Mix')

The Quarterly rebalancing with 10% threshold produces an attractive portfolio with top performance and low turnover.

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

Extending Gold time series

September 11, 2012 2 comments

While back-testing trading strategies I want all assets to have long history. Unfortunately, sometimes there is no tradeable stock or ETF with sufficient history. For example, I might use GLD as a proxy for Gold allocation, but GLD is only began trading in November of 2004.

We can extend the GLD’s historical returns with its benchmark index – London Gold afternoon fixing prices. Let’s first download and plot side by side the GLD and London Gold afternoon fixing prices from http://www.kitco.com/gold.londonfix.html the convenient CSV download is provided by Wikiposit.

###############################################################################
# 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
	#****************************************************************** 
	load.packages('quantmod')	
	GLD = getSymbols('GLD', src = 'yahoo', from = '1970-01-01', auto.assign = F)			
		GLD = adjustOHLC(GLD, use.Adjusted=T)
		
	# get Gold.PM - London Gold afternoon fixing prices
	temp = read.csv('http://wikiposit.org/w?action=dl&dltypes=comma%20separated&sp=daily&uid=KITCO',skip=4,header=TRUE, stringsAsFactors=F)
	Gold.PM = make.xts(as.double(temp$Gold.PM) / 10, as.Date(temp$Date, '%d-%b-%Y'))
		Gold.PM = Gold.PM[ !is.na(Gold.PM), ]

	# merge GLD and Gold.PM
	data <- new.env()
		data$GLD = GLD
		data$Gold.PM = Gold.PM
	bt.prep(data, align='remove.na')
	
	#*****************************************************************
	# Plot GLD and Gold.PM
	#****************************************************************** 
	layout(1:2)
	plota(data$GLD, type='l', col='black', plotX=F)	
		plota.lines(data$Gold.PM, col='blue')	
	plota.legend('GLD,Gold.PM', 'black,blue', list(data$GLD, data$Gold.PM))
		
	# plot GLD and London afternoon gold fix spread
	spread = 100 * (Cl(data$GLD) - data$Gold.PM) / data$Gold.PM
	plota(spread , type='l', col='black')	
		plota.legend('GLD vs Gold.PM % spread', 'black', spread)

There is a divergence between GLD and its benchmark index. Please read “Gold’s ‘Paper’ Price” article by B. Zigler for a detailed discussion and explanation for this divergence.

Next let’s extend the GLD’s time series with its benchmark ( I created a helper function extend.GLD() function in data.r at github to simplify the process ) and use it for a simple equal weight strategy

	#*****************************************************************
	# Create simple equal weight back-test
	#****************************************************************** 
	tickers = spl('GLD,TLT')

	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)
		
		# extend GLD with Gold.PM - London Gold afternoon fixing prices
		data$GLD = extend.GLD(data$GLD)
		
	bt.prep(data, align='remove.na')
 
    #*****************************************************************
    # Code Strategies
    #******************************************************************
    prices = data$prices      
    n = ncol(prices)
  
    # find period ends
    period.ends = endpoints(prices, 'months')
        period.ends = period.ends[period.ends > 0]
        
    models = list()
   
    #*****************************************************************
    # Equal Weight
    #******************************************************************
    data$weight[] = NA
        data$weight[period.ends,] = ntop(prices[period.ends,], n)   
    models$equal.weight = bt.run.share(data, clean.signal=F)
    
    #*****************************************************************
    # Create Report
    #******************************************************************       
    plotbt.custom.report.part1(models)       
    plotbt.custom.report.part2(models)       

By extending GLD’s time series, the back-test now goes back to the June 2002, the TLT’s first trading date.

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

Categories: Backtesting, R
Follow

Get every new post delivered to your Inbox.

Join 246 other followers