[R-SIG-Finance] Help replicating a paper

[Ilya Kipnis]

I'm not sure this question gets posed often on this list, but here goes...

I'm trying to emulate the results from SSRN id2450017 aka Momentum,
Markowitz, and Smart Beta by W.J. Keller. My issue is that for a weighting
algorithm, it often only invests in one assets, no assets, or very few
assets. Furthermore, while using the 9 sector spiders, my results are
pretty terrible, and certainly nothing worthy of publishing in a paper, so
I know there's a mistake on my end somewhere, but I'm not quite sure where.

Here's the link to the paper:

http://papers.ssrn.com/sol3/papers.cfm?abstract_id=2450017

And my implementation of the MAA algorithm (special case at the moment,
will generalize later). I'm hoping my comments will help make the code
easier to follow.

*****************************

require(quantmod)
require(PerformanceAnalytics)

#Use 9 sector spiders
suffix <- c("Y", "P", "E", "F", "V", "I", "B", "K", "U")
symbols <- paste0("XL", suffix)
getSymbols(symbols, from="1990-01-01")
prices <- list()
for(i in 1:length(symbols)) {
  prices[[i]] <- Ad(get(symbols[i]))
}
prices <- do.call(cbind, prices)

rets <- Return.calculate(prices)
rets <- rets[-1,] #remove first row of NAs

weights <- list()
dates <- list()
monthlyEndpoints <- endpoints(rets, on = "months")
for(i in 2:(length(monthlyEndpoints) - 4)) {

  #subset returns and define mktRets as average daily returns
  returnsData <- rets[monthlyEndpoints[i]:monthlyEndpoints[i+4],]
  mktRets <- xts(rowMeans(returnsData), order.by=index(returnsData))

  #raw computed values -- annualized variance, cumulative returns
  #correlation to average daily returns, and annualized variance of
  #average daily returns
  variances <- as.numeric(StdDev.annualized(returnsData)^2)
  cumRets <- as.numeric(Return.cumulative(returnsData))
  cors <- as.numeric(cor(returnsData, mktRets))
  mktVar <- as.numeric(StdDev.annualized(mktRets) ^ 2)

  #shrink asset variances to the average by 50%
  #shrink cumulative returns to the average by 50%
  #shrink correlations to the average by 50%
  #shrink market variance to zero by 50%
  #see Keller (June 2014 Table 1 pg. 8)
  shrunkVariances <- (variances + mean(variances)) / 2
  shrunkCumRets <- (cumRets + mean(cumRets)) / 2
  shrunkCors <- (cors+mean(cors))/2
  shrunkMktVar <- mktVar/2

  #values computed from shrunken values
  #define beta as shrunken variances * shrunken correlations / shrunken
market variance
  #define idiosyncratic variances as shrunken variances minus shrunken
market variance * beta
  #see Keller (June 2014 under A.4)
  betas <- shrunkVariances*shrunkCors/shrunkMktVar
  idiosyncraticVars <- shrunkVariances-shrunkMktVar*betas

  #compute asset treynors
  assetTreynors <- shrunkCumRets/betas

  #set initial threshold to zero
  longs <- assetTreynors > 0

  converge <- FALSE
  while(!converge) {
    old <- longs

    #implement equation A.6 from Keller (June 2014)
    newTreynorThreshNum <- shrunkMktVar * sum(old * shrunkCumRets * betas /
idiosyncraticVars)
    newTreynorThreshDenom <- 1 + shrunkMktVar * sum(old * betas * betas /
idiosyncraticVars)
    newTreynorThresh <- newTreynorThreshNum/newTreynorThreshDenom

    #assets remaining are those that have a higher initial treynor ratio
    longs <- assetTreynors > newTreynorThresh

    #if our assets are the same as they were before the computation of the
new threshold, end
    if(sum(old - longs)==0) {
      converge <- TRUE
    }
  }

  #after convergence, allocate weights according to w_i~(1-t/t_i)*r_i/s_i
and normalize
  longs <-
longs*(1-newTreynorThresh/assetTreynors)*shrunkCumRets/idiosyncraticVars
  longs <- longs/sum(longs)

  weights[[i]] <-longs
  dates[[i]] <- as.character(index(returnsData)[nrow(returnsData)])
}

#set up weights xts
weights <- do.call(rbind, weights)
dates <- do.call(rbind, dates)
weights <- xts(weights, order.by=as.Date(dates))
weights[is.na(weights)] <- 0
colnames(weights) <- colnames(rets)

#for dates without investment, invest in risk-free
rets$riskFree <- 0
weights$riskFree <- rep(1, nrow(weights))-rowSums(weights)

#results are pretty terrible
test <- Return.rebalancing(R = rets, weights = weights)
charts.PerformanceSummary(test)
plot(log(cumprod(1+test)))
Return.annualized(test) #less than 5%
maxDrawdown(test) #greater than 30%

#Often invests in one security
head(round(weights, 3), 20)

*******************************
Thank you for any help ahead of time.

-Ilya

	[[alternative HTML version deleted]]