mcsim gist for rapid prototypive collaborative narrative R

mcsim_gist.R

#' Monte Carlo simulate strategy results
#'
#' Return bands of returns based on Monte Carlo simulations of back-test results
#' @param Portfolio string identifier of portfolio name
#' @param Account string identifier of account name
#' @param n number of simulations, default = 1000
#' @param \dots any other passthrough parameters
#' @param l block length, default = 20
#' @param use determines whether to use 'daily' or 'txn' PL, default = "equity" ie. daily
#' @param gap numeric number of periods from start of series to start on, to eliminate leading NA's
#' @return a plot.xts object of simulated return streams
#' @return a list containing:
#' \itemize{
#'   \item{\code{replicates}:}{an xts object containing all the resampled time series replicates}
#'   \item{\code{dailypl}:}{an xts object containing daily P&L from the original backtest}
#'   \item{\code{initeq}:}{a numeric variable containing the initEq of the portfolio, for starting portfolio value}
#'   \item{\code{num}:}{a numeric variable reportin the number of replicaes in the simulation}
#'   \item{\code{length}:}{a numeric variable reporting the block length used in the simulation, if any}
#'   \item{\code{ddvec}:}{a numeric vector of drawdowns for each replicate series}
#'   \item{\code{w}:}{a string containing information on whether the simulation is with or without replacement}
#'   \item{\code{call}:}{an object fo type \code{call} that contains the \code{mcsim} call}
#' }
#'
#' Note that this object and its slots may change in the future.
#' Slots \cod{replicates},\code{dailypl},\code{initeq}, and \code{call} are likely
#' to exist in all future versions of this function, but other slots may be added
#' and removed as \code{S3method}s are developed.
#'
#' @note
#' Requires boot package
#' @importFrom boot tsboot
#' @importFrom boot boot.array
#' @author Jasen Mackie, Brian G. Peterson
#' @seealso \code{\link{boot}}
#' @export
mcsim <- function(  Portfolio
                    , Account
                    , n = 1000
                    , replacement = TRUE
                    , ...
                    , l = 1
                    , use=c('equity','txns')
                    , gap = 1
){
  
  use=use[1] #take the first value if the user didn't specify
  switch (use,
          Eq =, eq =, Equity =, equity =, cumPL = {
            dailyPL <- dailyEqPL(Portfolio, incl.total = TRUE)
            dailyPL <- dailyPL[gap:nrow(dailyPL), ncol(dailyPL)]
          },
          Txns =, txns =, Trades =, trades = {
            dailyPL <- dailyTxnPL(Portfolio,  incl.total = TRUE)
            dailyPL <- dailyPL[gap:nrow(dailyPL), ncol(dailyPL)]
          }
  )
  
  p <- getPortfolio(Portfolio)
  a <- getAccount(Account)
  initEq <- attributes(a)$initEq
  use=c('equity','txns')
  tmp <- NULL
  k <- NULL
  EndEqdf <- data.frame(dailyPL)
  
  if(isTRUE(replacement)) {
    #tsb <- tsboot(ROC(initEq + cumsum(dailyPL)), maxDrawdown, invert=FALSE, n, l, sim = "fixed", ...)
    tsb <- tsboot(dailyPL, function(x) { -max(cummax(cumsum(x))-cumsum(x)) }, n, l, sim = "fixed", ...)
    inds <- t(boot.array(tsb))
    #k <- NULL
    tsbootARR <- NULL
    tsbootxts <- NULL
    EndEqdf[is.na(EndEqdf)] <- 0
    for(k in 1:ncol(inds)){
      tmp <- cbind(tmp, EndEqdf[inds[,k],])
    }
    tsbootARR <- apply(tmp, 2, function(x) cumsum(x))
    which(is.na(tsbootARR))
    tsbootxts <- xts(tsbootARR, index(dailyPL))
    dd <- tsb$t
    withorwithout <- "WITH REPLACEMENT"
  } else if(!isTRUE(replacement)) {
    tsbootxts <- foreach (k=1:n, .combine=cbind.xts ) %dopar% {
      if(isTRUE(l>1)) {
#         # do a block resample, without replacement
#         # first, resample the index
#         x <- 1:length(dailyPL)
#         # now sample blocks
#         # TODO: this method creates variable length blocks centered on 'l'
#         #       we could also do fixed blocks with seq_along() to create 's'
#         s <- sort(sample(x=x[2:length(x)-1],size = l-1,replace = FALSE))
#         blocks<-split(x, findInterval(x,s))
        x <- 1:length(dailyPL) # uneven total length
        #l<-10
        blocks <- seq(sample.int(1:l,1),length(x),by=l)
        # now reassemble the target index order
        idx <- unlist(blocks[sample(names(blocks),size = length(blocks),replace = FALSE)]) ; names(idx)<-NULL
        tmp <- as.vector(dailyPL)[idx]
      } else {
        tmp <- sample(as.vector(dailyPL), replace = replacement)
      }
      tsbootARR <- cumsum(tmp)
      tsbootxts <- xts(tsbootARR, index(dailyPL))
    }
    dd <- apply(tsbootxts, 2, function(x) { -max(cummax(x)-(x)) } )
    withorwithout <- "WITHOUT REPLACEMENT"
  }
  
  ret <- list(replicates=tsbootxts, dailypl=dailyPL, initeq=initEq, num=n, length=l, ddvec=dd,
              w=withorwithout, call=match.call())
  class(ret) <- "mcsim"
  ret
}



#' plot method for objects of type \code{mcsim}
#'
#' @param x object of type 'mcsim' to plot
#' @param y not used, to match generic signature, may hold overlay daya in the future
#' @param \dots any other passthrough parameters
#' @param normalize TRUE/FALSE whether to normalize the plot by initEq, default TRUE
#' @author Jasen Mackie, Brian G. Peterson
#'
#' @export
plot.mcsim <- function(x, y, ..., normalize=TRUE) {
  ret <- x
  if(isTRUE(normalize) && ret$initeq>1 ){
    div <- ret$initeq
  } else {
    div <- 1
  }
  p   <- plot.xts(ret$replicates/div
                  , col = "lightgray"
                  , main = paste0("Sample Backtest ",ret$w)
                  , grid.ticks.on = 'years'
  )
  p   <- lines(cumsum(ret$dailypl/div), col = "red")
  p
}

#' hist method for objects of type \code{mcsim}
#'
#' @param x object of type mcsim to plot
#' @param \dots any other passthrough parameters
#' @author Jasen Mackie, Brian G. Peterson
#'
#' @export
hist.mcsim <- function(x, ..., normalize=TRUE) {
  ret <- x
  if(isTRUE(normalize) && ret$initeq>1 ){
    # running maxDD vector - note that 1st day is zero...need to fix
    idxarr <- -(cummax(cumsum(ret$dailypl))-cumsum(ret$dailypl))
    # maxDD from idxarr
    maxDD <- min(idxarr)
    # find index of maxDD
    idxmaxdd <- which(maxDD == idxarr)[1] # subset 1st element if vector length > 1
    # add max cum PL with initeq to get denom to use for normalizing maxDD
    divhist <- cummax(cumsum((ret$dailypl)))[idxmaxdd]+ret$initeq
    
    # now calculate div for sampled maxDDs
    idxarrtsb <- apply(ret$replicates, 2, function(x) { -(cummax(x)-(x)) } )
    maxDDtsb <- apply(idxarrtsb, 2, function(x) {min(x)} )
    idxmaxddtsb <- NULL
    for (i in 1:ret$num) {
      idxmaxddtsb[i] <- which(maxDDtsb[i] == idxarrtsb[,i])[1]
    }
    div <- NULL
    for (i in 1:ret$num) {
      div[i] <- cummax(ret$replicates[,i])[idxmaxddtsb[i]] + ret$initeq
    }
    
    # and finally calculate div for sample median
    divmed <- NULL
    divmed <- ifelse(ret$num%%2!=0, div[which(ret$ddvec==median(ret$ddvec))], div[which(ret$ddvec[-1]==median(ret$ddvec[-1]))])

    
  } else {
    divhist <- 1
    div <- 1
    divmed <- 1
  }
  xname <- paste(ret$num, "replicates with block length", ret$length)
  h <- hist(ret$ddvec/div, main = paste("Drawdown distribution", ret$w, "of" , xname), breaks="FD"
            , xlab = "Max Drawdown", ylab = "Density"
            , col = "lightgray", border = "white", freq=FALSE
  )
  h
  box(col = "darkgray")
  
  b = -max(cummax(cumsum(ret$dailypl))-cumsum(ret$dailypl))/divhist
  #b = maxDrawdown(ROC(ret$initeq + cumsum(ret$dailypl)), invert = FALSE)
  abline(v = b, col = "darkgray", lty = 2)
  b.label = ("Backtest Max Drawdown")
  h = rep(0.2 * par("usr")[3] + 1 * par("usr")[4], length(b))
  text(b, h, b.label, offset = 0.2, pos = 2, cex = 0.8, srt = 90)
  
  abline(v=median(na.omit(ret$ddvec))/divmed, col = "darkgray", lty = 2)
  c.label = ("Sample Median Max Drawdown")
  text(median(na.omit(ret$ddvec))/divmed, h, c.label, offset = 0.2, pos = 2, cex = 0.8, srt = 90)
}

#' quatile method for objects of type \code{mcsim}
#'
#' @param x object of type 'mcsim' to produce replicate quantiles
#' @param \dots any other passthrough parameters
#' @author Jasen Mackie, Brian G. Peterson
#'
#' @export
quantile.mcsim <- function(x, ..., normalize=TRUE) {
  ret <- x
  q   <- quantile(ret$replicates)
  q
}

mysim <- mcsim("longtrend", "longtrend", n=1000, replacement=FALSE, l=10, gap=10)
plot.mcsim(mysim, normalize = FALSE)
hist.mcsim(mysim, normalize = TRUE)
hist.mcsim(mysim, normalize = FALSE)
quantile.mcsim(mysim)

mysim <- mcsim("longtrend", "longtrend", n=1001, replacement=TRUE, l=30, gap=1)
plot.mcsim(mysim, normalize = FALSE)
hist.mcsim(mysim, normalize = TRUE)
hist.mcsim(mysim, normalize = FALSE)
quantile.mcsim(mysim)
# 
# tailmysim <- mcsim("bbands", "bbands", n=10, replacement=TRUE, l=10, gap=10)
# plot.mcsim(tailmysim)
# hist.mcsim(tailmysim)
# 
# headmysim <- mcsim("longtrend", "longtrend", n=10, replacement=TRUE, l=10, gap=10)
# plot.mcsim(headmysim)
# hist.mcsim(headmysim)

###############################################################################
# Blotter: Tools for transaction-oriented trading systems development
# for R (see http://r-project.org/)
# Copyright (c) 2008-2016 Peter Carl and Brian G. Peterson
#
# This library is distributed under the terms of the GNU Public License (GPL)
# for full details see the file COPYING
#
# $Id$
#
###############################################################################