Archive
Fast Threshold Clustering Algorithm (FTCA) test
Today I want to share the test and implementation for the Fast Threshold Clustering Algorithm (FTCA) created by David Varadi. This implementation was developed and contributed by Pierre Chretien, I only made minor updates.
Let’s first replicate the results from the Fast Threshold Clustering Algorithm (FTCA) post:
###############################################################################
# Load Systematic Investor Toolbox (SIT)
# https://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('XLY,XLP,XLE,XLF,XLV,XLI,XLB,XLK,XLU')
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='keep.all')
#*****************************************************************
# Helper function to compute portfolio allocation additional stats
#******************************************************************
portfolio.allocation.custom.stats.clusters <- function(x,ia) {
return(list(
clusters.FTCA = cluster.group.FTCA(0.5)(ia)
))
}
#*****************************************************************
# Find clusters
#******************************************************************
periodicity = 'months'
lookback.len = 252
obj = portfolio.allocation.helper(data$prices,
periodicity = periodicity, lookback.len = lookback.len,
min.risk.fns = list(EW=equal.weight.portfolio),
custom.stats.fn = portfolio.allocation.custom.stats.clusters
)
clusters = obj$clusters.FTCA$EW
clusters['2012:05::']
The clusters are stable and match David’s results
XLB XLE XLF XLI XLK XLP XLU XLV XLY
2012-05-31 1 1 1 1 1 1 1 1 1
2012-06-29 1 1 1 1 1 1 1 1 1
2012-07-31 1 1 1 1 1 1 1 1 1
2012-08-31 1 1 1 1 1 1 1 1 1
2012-09-28 1 1 1 1 1 1 1 1 1
2012-10-31 1 1 1 1 1 1 1 1 1
2012-11-30 2 2 2 2 2 2 1 2 2
2012-12-31 2 2 2 2 2 2 1 2 2
2013-01-31 2 2 2 2 2 2 1 2 2
2013-02-28 1 1 1 1 1 1 1 1 1
2013-03-28 1 1 1 1 1 1 1 1 1
2013-04-30 1 1 1 1 1 1 1 1 1
2013-05-31 1 1 1 1 1 1 1 1 1
2013-06-28 1 1 1 1 1 1 1 1 1
2013-07-31 1 1 1 1 1 1 1 1 1
2013-08-30 1 1 1 1 1 1 1 1 1
2013-09-30 1 1 1 1 1 1 1 1 1
2013-10-31 1 1 1 1 1 1 1 1 1
2013-11-26 1 1 1 1 1 1 1 1 1
Next let’s compare the Cluster Portfolio Allocation Algorithm using K-means and FTCA:
#***************************************************************** # Code Strategies #****************************************************************** obj = portfolio.allocation.helper(data$prices, periodicity = periodicity, lookback.len = lookback.len, min.risk.fns = list( C.EW.kmeans = distribute.weights(equal.weight.portfolio, cluster.group.kmeans.90), C.EW.FTCA = distribute.weights(equal.weight.portfolio, cluster.group.FTCA(0.5)) ) ) models = create.strategies(obj, data)$models #***************************************************************** # Create Report #****************************************************************** barplot.with.labels(sapply(models, compute.turnover, data), 'Average Annual Portfolio Turnover')
Both clustering algorithms produced very similar results. One noticeable difference is turnover. Since the Fast Threshold Clustering Algorithm (FTCA) produced more stable groups, it had smaller turnover.
The full source code and example for the cluster.group.FTCA() function is available in strategy.r at github.
getSymbols Extra
The getSymbols function from the quantmod package is an easy and convenient way to bring historical stock prices into your R environment. You need to specify the list of tickers, the source of historical prices and dates. For example following commands will download historical stock prices from yahoo finance for ‘RWX’, ‘VNQ’, ‘VGSIX’ symbols:
data <- new.env()
getSymbols.extra(c('RWX','VNQ','VGSIX'), src = 'yahoo', from = '1980-01-01', env = data, auto.assign = T)
Now, the data environment contains the historical stock prices. For example, create plot for RWX using the code below:
plot(data$RWX)
Sometimes, I find that getSymbols functionality can be extended. For example, it would be nice to rename the series. I.e. RWX is a real estate ETF, so we could ask getSymbols function to get RWX, but call output Real.Estate. Another useful feature would be ability to specify how to extend the data. I.e. RWX is only started trading in 2007, it would be convenient to extend RWX time series with VNQ and VGSIX.
I created the getSymbols.extra function to address these features. The getSymbols.extra function allows you to specify tickers in the following format:
- RWX : the original functionality. i.e get historical stock prices for RWX and it can be accessed with data$RWX
- Real.Estate = RWX : get historical stock prices for RWX, rename it to Real.Estate, and it can be accessed with data$Real.Estate
- RWX + VNQ + VGSIX : get historical stock prices for RWX, VNQ, VGSIX and extend RWX with VNQ, next extend it with VGSIX, and it can be accessed with data$RWX
- Real.Estate = RWX + VNQ + VGSIX : mix and match above functionality
Let’s look at the example:
###############################################################################
# Load Systematic Investor Toolbox (SIT)
# https://systematicinvestor.wordpress.com/systematic-investor-toolbox/
###############################################################################
setInternet2(TRUE)
con = gzcon(url('http://www.systematicportfolio.com/sit.gz', 'rb'))
source(con)
close(con)
tickers = spl('REIT=RWX, RWX+VNQ, REIT.LONG=RWX+VNQ+VGSIX')
data <- new.env()
getSymbols.extra(tickers, src = 'yahoo', from = '1980-01-01', env = data, auto.assign = T)
bt.start.dates(data)
data$symbolnames = spl('REIT.LONG,RWX,REIT')
for(i in data$symbolnames) data[[i]] = adjustOHLC(data[[i]], use.Adjusted=T)
bt.prep(data, align='keep.all')
plota.matplot(data$prices)
The REIT.LONG is RWX extended using VNQ and next using VGSIX.
There is also a possibility to provide custom data into the getSymbols.extra function. I.e. the data which is not available through the getSymbols function.
For example, to extend GLD time series we might used historical Gold prices from Bundes Bank:
# Use extrenal data
raw.data <- new.env()
raw.data$GOLD = bundes.bank.data.gold()
tickers = spl('GLD, GLD.LONG=GLD+GOLD')
data <- new.env()
getSymbols.extra(tickers, src = 'yahoo', from = '1980-01-01', env = data, raw.data = raw.data, auto.assign = T)
bt.start.dates(data)
data$symbolnames = spl('GLD.LONG,GLD')
for(i in data$symbolnames) data[[i]] = adjustOHLC(data[[i]], use.Adjusted=T)
bt.prep(data, align='keep.all')
plota.matplot(data$prices)
The main goal of the getSymbols.extra function is to simplify your data acquisition/preparation step. Please let me know if you run into any problems or have suggestions.
To view the complete source code for this example, please have a look at the getSymbols.extra.test() function in utils.r at github.
Running Back-tests in parallel
Once you start experimenting with many different asset allocation algorithms, the computation time of running the back-tests can be substantial. One simple way to solve the computation time problem is to run the back-tests in parallel. I.e. if the asset allocation algorithm does not use the prior period holdings to make decision about current allocation, we can run many periods in parallel.
In the Update for Backtesting Asset Allocation Portfolios post, I show cased the portfolio.allocation.helper() function in strategy.r at github. The portfolio.allocation.helper() function is a user-friendly interface to evaluate multiple asset allocation algorithms over given asset universe in a sequential fashion.
Following is a sample code from the Update for Backtesting Asset Allocation Portfolios post:
#*****************************************************************
# Code Strategies
#******************************************************************
obj = portfolio.allocation.helper(data$prices,
periodicity = 'months', lookback.len = 60,
min.risk.fns = list(
EW=equal.weight.portfolio,
RP=risk.parity.portfolio,
MD=max.div.portfolio,
MV=min.var.portfolio,
MVE=min.var.excel.portfolio,
MV2=min.var2.portfolio,
MC=min.corr.portfolio,
MCE=min.corr.excel.portfolio,
MC2=min.corr2.portfolio,
MS=max.sharpe.portfolio(),
ERC = equal.risk.contribution.portfolio,
# target retunr / risk
TRET.12 = target.return.portfolio(12/100),
TRISK.10 = target.risk.portfolio(10/100),
# rso
RSO.RP.5 = rso.portfolio(risk.parity.portfolio, 5, 500),
# others
MMaxLoss = min.maxloss.portfolio,
MMad = min.mad.portfolio,
MCVaR = min.cvar.portfolio,
MCDaR = min.cdar.portfolio,
MMadDown = min.mad.downside.portfolio,
MRiskDown = min.risk.downside.portfolio,
MCorCov = min.cor.insteadof.cov.portfolio
)
)
To run the same strategies in parallel, I created the portfolio.allocation.helper.parallel() function in strategy.r at github. There is one extra input that you need to specify: cores – number of CPU processors used for computations.
For example, the code below will use 2 CPU processors to run back-test computations. It will run faster than the portfolio.allocation.helper() function.
#*****************************************************************
# Code Strategies
#******************************************************************
obj = portfolio.allocation.helper.parallel(cores = 2, data$prices,
periodicity = 'months', lookback.len = 60,
min.risk.fns = list(
EW=equal.weight.portfolio,
RP=risk.parity.portfolio,
MD=max.div.portfolio,
MV=min.var.portfolio,
MVE=min.var.excel.portfolio,
MV2=min.var2.portfolio,
MC=min.corr.portfolio,
MCE=min.corr.excel.portfolio,
MC2=min.corr2.portfolio,
MS=max.sharpe.portfolio(),
ERC = equal.risk.contribution.portfolio,
# target retunr / risk
TRET.12 = target.return.portfolio(12/100),
TRISK.10 = target.risk.portfolio(10/100),
# rso
RSO.RP.5 = rso.portfolio(risk.parity.portfolio, 5, 500),
# others
MMaxLoss = min.maxloss.portfolio,
MMad = min.mad.portfolio,
MCVaR = min.cvar.portfolio,
MCDaR = min.cdar.portfolio,
MMadDown = min.mad.downside.portfolio,
MRiskDown = min.risk.downside.portfolio,
MCorCov = min.cor.insteadof.cov.portfolio
)
)
Hopefully, I did not ruin your prolong lunch plans 🙂
Commissions
Today, I want to explain the commission’s functionality build in to Systematic Investor Toolbox(SIT) “share” back-test.
At each re-balance time the capital is allocated given the weight such that
share = weight * capital / price cash = capital - share * price
For example, if weight is 100% (i.e. fully invested) and capital = $100 and price = $10 then
share = 10 shares cash = $0
The period return is equal to
return = [share * price + cash] / [share * price.yesterday + cash] - 1
The total return is equal to
total.return = [1 + return.0] * [1 + return.1] * ... * [1 + return.n] - 1
The period returns constructed this way let me construct portfolio returns without using loops
I.e. if share, price, cash are matrices, then
portfolio.return = rowSums((share * price + cash) / (share * mlag(price) + cash) - 1)
and
equity = cumprod(1 + portfolio.return)
To introduce commissions into above framework, I had to make following assumptions.
Let’ assume that P0 and P1 are stock prices and com is commission that is very small relative to stock price then
[P0 - com] / P0 is close (equal) to P0 / [P0 + com] and [P0 - com] * P1 is close (equal) to P0 * [P1 - com]
Given commissions, the period return formula used in SIT is equal to
return = [share * price + cash - commission] / [share * price.yesterday + cash] - 1
Now let’s look at the example of trade with commissions:
Let's say we are fully invested (i.e. cash = 0 and capital = share * P0) opening trade cost = share * P0 + com closing trade cost = share * P1 - com return = [closing trade cost] / [opening trade cost] - 1 = [share * P1 - com] / [share * P0 + com] - 1
In SIT, these computations are equivalent to
([capital - com] / [capital]) * ([share * P1 + cash - com] / [share * P0 + cash]) - 1 < given that cash = 0 and capital = share * P0 > = ([share * P0 - com] / [share * P0]) * ([share * P1 - com] / [share * P0]) - 1 < given [P0 - com] / P0 ~ P0 / [P0 + com] > = ([share * P0] / [share * P0 + com] / ) * ([share * P1 - com] / [share * P0]) - 1 = [share * P1 - com] / [share * P0 + com] - 1
Hence, as long as commissions are small relative to whole trade the returns produced with SIT will be very close to the true returns.
SIT currently supports following 3 types of commissions
- cents / share commission
commission = abs(share - mlag(share)) * commission$cps
- fixed commission
commission = sign(abs(share - mlag(share))) * commission$fixed
- percentage commission
commission = price * abs(share - mlag(share)) * commission$percentage
You can mix and match these commission methods in any way,
the total.commission = cps.commission + fixed.commission + percentage.commission
and period return is equal to
return = (share * price + cash - total.commission) / (share * mlag(price) + cash) - 1
Next let’s see the impact of different type of commissions. There two ways to specify the commissions.
-
commission = 0.1
will be translated in
commission = list(cps = 0.1, fixed = 0.0, percentage = 0.0)
-
commission = list(cps = 0.0, fixed = 0.0, percentage = 0.0)
###############################################################################
# Load Systematic Investor Toolbox (SIT)
# https://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('EEM')
data <- new.env()
getSymbols(tickers, src = 'yahoo', from = '1970-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='2013:08::2013:09')
#*****************************************************************
# Code Strategies
#******************************************************************
buy.date = '2013:08:14'
sell.date = '2013:08:15'
day.after.sell.date = '2013:08:16'
capital = 100000
prices = data$prices
share = as.double(capital / prices[buy.date])
# helper function to compute trade return
comp.ret <- function(sell.trade.cost, buy.trade.cost) { round(100 * (as.double(sell.trade.cost) / as.double(buy.trade.cost) - 1), 2) }
#*****************************************************************
# Zero commission
#******************************************************************
data$weight[] = NA
data$weight[buy.date] = 1
data$weight[sell.date] = 0
commission = 0.0
model = bt.run.share(data, commission = commission, capital = capital, silent = T)
comp.ret( share * prices[sell.date], share * prices[buy.date] )
comp.ret( model$equity[day.after.sell.date], model$equity[buy.date] )
#*****************************************************************
# 10c cps commission
# cents / share commission
# trade cost = abs(share - mlag(share)) * commission$cps
#******************************************************************
data$weight[] = NA
data$weight[buy.date] = 1
data$weight[sell.date] = 0
commission = 0.1
model = bt.run.share(data, commission = commission, capital = capital, silent = T)
comp.ret( share * (prices[sell.date] - commission), share * (prices[buy.date] + commission) )
comp.ret( model$equity[day.after.sell.date], model$equity[buy.date] )
#*****************************************************************
# $5 fixed commission
# fixed commission per trade to more effectively to penalize for turnover
# trade cost = sign(abs(share - mlag(share))) * commission$fixed
#******************************************************************
data$weight[] = NA
data$weight[buy.date] = 1
data$weight[sell.date] = 0
commission = list(cps = 0.0, fixed = 5.0, percentage = 0.0)
model = bt.run.share(data, commission = commission, capital = capital, silent = T)
comp.ret( share * prices[sell.date] - commission$fixed, share * prices[buy.date] + commission$fixed )
comp.ret( model$equity[day.after.sell.date], model$equity[buy.date] )
#*****************************************************************
# % commission
# percentage commission
# trade cost = price * abs(share - mlag(share)) * commission$percentage
#******************************************************************
data$weight[] = NA
data$weight[buy.date] = 1
data$weight[sell.date] = 0
commission = list(cps = 0.0, fixed = 0.0, percentage = 1/100)
model = bt.run.share(data, commission = commission, capital = capital, silent = T)
comp.ret( share * prices[sell.date] * (1 - commission$percentage), share * prices[buy.date] * (1 + commission$percentage) )
comp.ret( model$equity[day.after.sell.date], model$equity[buy.date] )
To view the complete source code for this example, please have a look at the bt.commission.test() function in bt.test.r at github.
Weekend Reading: Market Neutral
I recently came across a very interesting idea at the The Problem with Market Neutral (and an Answer) post by Mebane Faber. Today I want to show how you can test such strategy using the Systematic Investor Toolbox:
###############################################################################
# Load Systematic Investor Toolbox (SIT)
# https://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')
data = new.env()
# load historical market returns
temp = get.fama.french.data('F-F_Research_Data_Factors', periodicity = '',download = T, clean = T)
ret = temp[[1]]$Mkt.RF + temp[[1]]$RF
price = bt.apply.matrix(ret / 100, function(x) cumprod(1 + x))
data$SPY = make.stock.xts( price )
# load historical momentum returns
temp = get.fama.french.data('10_Portfolios_Prior_12_2', periodicity = '',download = T, clean = T)
ret = temp[[1]]
price = bt.apply.matrix(ret / 100, function(x) cumprod(1 + x))
data$HI.MO = make.stock.xts( price$High )
data$LO.MO = make.stock.xts( price$Low )
# align dates
bt.prep(data, align='remove.na')
#*****************************************************************
# Code Strategies
#*****************************************************************
models = list()
data$weight[] = NA
data$weight$SPY[] = 1
models$SPY = bt.run.share(data, clean.signal=T)
data$weight[] = NA
data$weight$HI.MO[] = 1
models$HI.MO = bt.run.share(data, clean.signal=T)
data$weight[] = NA
data$weight$LO.MO[] = 1
models$LO.MO = bt.run.share(data, clean.signal=T)
data$weight[] = NA
data$weight$HI.MO[] = 1
data$weight$LO.MO[] = -1
models$MKT.NEUTRAL = bt.run.share(data, clean.signal=F)
#*****************************************************************
# Modified MN
# The modified strategy below starts 100% market neutral, and depending on the drawdown bucket
# will reduce the shorts all the way to zero once the market has declined by 50%
# (in 20% steps for every 10% decline in stocks)
#*****************************************************************
market.drawdown = -100 * compute.drawdown(data$prices$SPY)
market.drawdown.10.step = 10 * floor(market.drawdown / 10)
short.allocation = 100 - market.drawdown.10.step * 2
short.allocation[ short.allocation < 0 ] = 0
data$weight[] = NA
data$weight$HI.MO[] = 1
data$weight$LO.MO[] = -1 * short.allocation / 100
models$Modified.MN = bt.run.share(data, clean.signal=F)
#*****************************************************************
# Create Report
#*****************************************************************
strategy.performance.snapshoot(models, T)
Mebane thank you very much for sharing this great observation and great strategy that works! I would encourage readers to experiment with idea and share their findings.
If you want to concentrate on the long side, one idea that comes to mind is to start not fully invested say at 90% allocation, and once the market hits say 20% draw-down to invest 100% in expectation of quick recovery.
To view the complete source code for this example, please have a look at the bt.mebanefaber.modified.mn.test() function in bt.test.r at github.
Systematic Investor 