Backtesting Part 5

Backtesting Part 5.R

rm(list = ls(all = TRUE))
memory.limit(size = 4095)

#Get Data
library(quantmod)
getSymbols('^GSPC',from='1900-01-01')
myStock <- Cl(GSPC)
bmkReturns <- dailyReturn(myStock, type = "arithmetic")

#Apply our strategy to tomorrows returns
#Today's close to tomorrow's close
myReturns <- Next(bmkReturns)
myReturns[nrow(myReturns)] <- 0 

#Calculate days since the n-day high
library(compiler)
daysSinceHigh <- cmpfun(function(x, n){
  apply(embed(x, n), 1, which.max)-1
})

#Setup parallel cluster
library(doParallel)
cl <- makeCluster(detectCores(), type='PSOCK')
registerDoParallel(cl)

#Calculate days since each n-day high for the stock
highs <- seq(5,500,by=5)
highMatrix <- matrix(data=NA,nrow=length(myStock),ncol=length(highs))
colnames(highMatrix) <- highs
system.time(
  highMatrix <- foreach(nHigh=highs) %dopar% {
    out <- daysSinceHigh(myStock,nHigh)
    out <- c(rep(NA,nHigh-1),out)
    return(out)
  })

highMatrix <- data.frame(do.call(cbind, highMatrix))
names(highMatrix) <- as.character(highs)

#Calculate Returns for various combinations of n-day highs and holding periods
#Thanks to G See for proposing the speedup
holds <- seq(5,500,by=5)
system.time(
  returnsList <- foreach(nHold=holds) %dopar% {
    out <- highMatrix
    out[out <= nHold] <- 1
    out[is.na(out) | out > nHold] <- 0
    sweep(out, MARGIN=1, myReturns, `*`)
  })
)
names(returnsList) <- holds

#Calculate Cumulative Returns for various scenarios
cumRet <- as.list(rep(NA,length(returnsList)))
i=1
for (returnMatrix in returnsList) {
  cumRet[[i]] <- apply(returnMatrix,MARGIN=2,function(x) prod(1 + x) - 1)
  i <- i+1
}
cumRet <- unlist(cumRet)
bmkRet <- prod(1 + bmkReturns) - 1
exRet <- cumRet-bmkRet

#Custom Color Ramp Function
range01 <- function(x)(x-min(x))/diff(range(x))
cRamp <- function(x){
  cols <- colorRamp(topo.colors(10))(range01(x))
  apply(cols, 1, function(xt)rgb(xt[1], xt[2], xt[3], maxColorValue=255))
}  

#Plot
Data <- data.frame(expand.grid(nHigh=highs,nHold=holds),exRet=exRet)
Data <- Data[Data$nHold<=Data$nHigh,]
plot(Data[,c(1,2)],col=cRamp(Data$exRet),pch=19,lwd=2)

#Stop parallel cluster
stopCluster(cl)