R/maxdd_portfolio.R

In rengelke/quantTraiding_trato: Trading Strategies Development Tool Kit

#' Calculate Optimized Portfolio with Drawdown Constraint
#'
#' @param R a xts object of asset returns
#' @param max_dd maximun drawdown to be calculated for the portfolio
#' @param soft_budget if TRUE, positions do not need to sum up to 100 percent
#' @param ... Additional arguments for Rglpk::Rglpk_solve_LP
#'
#' @return weights
#' @export
#' @importFrom Rglpk Rglpk_solve_LP
#'
#' @examples maxdd_portfolio(stock_returns, max_dd = 0.1, soft_budget = TRUE)
maxdd_portfolio <- function (R, max_dd = 0.1, soft_budget = FALSE, ...)  {

  if (is.null(dim(R))) {
    stop("Argument for 'R' must be rectangular.\\n")
  }
  if (any(is.na(R))) {
    stop("NA-values contained in object for 'R'.\\n")
  }
  if (max_dd <= 0 || max_dd >= 1) {
    stop("Argument for 'max_dd' must be in the interval (0, 1).\\n")
  }
  call <- match.call()
  RC <- as.matrix(cumsum(R))
  rownames(RC) <- NULL
  N <- ncol(RC)
  J <- nrow(RC)
  w <- rep(0, N)
  u <- rep(0, J)
  x <- c(w, u)
  obj <- c(as.numeric(RC[J, ]), rep(0, J))
  a1 <- cbind(diag(N), matrix(0, nrow = N, ncol = J))
  d1 <- rep(">=", N)
  b1 <- rep(0, N)
  a2 <- c(rep(1, N), rep(0, J))
  ifelse(soft_budget, d2 <- "<=", d2 <- "==")
  b2 <- 1
  a3 <- cbind(-1 * RC, diag(J))
  d3 <- rep("<=", J)
  b3 <- rep(max_dd, J)
  a4 <- a3
  d4 <- rep(">=", J)
  b4 <- rep(0, J)
  D1 <- -1 * diag(J)
  udiag <- embed(1:J, 2)[, c(2, 1)]
  D1[udiag] <- 1
  a5 <- cbind(matrix(0, ncol = N, nrow = J), D1)
  a5 <- a5[-J, ]
  d5 <- rep(">=", J - 1)
  b5 <- rep(0, J - 1)
  a6 <- c(rep(0, N), 1, rep(0, J - 1))
  d6 <- "=="
  b6 <- 0
  Amat <- rbind(a1, a2, a3, a4, a5, a6)
  Dvec <- c(d1, d2, d3, d4, d5, d6)
  Bvec <- c(b1, b2, b3, b4, b5, b6)
  opt <- Rglpk::Rglpk_solve_LP(obj = obj, mat = Amat, dir = Dvec,
                               rhs = Bvec, max = TRUE)
  if (opt$status != 0) {
    warning(paste("GLPK had exit status:", opt$status))
  }
  weights <- round(opt$solution[1:N], 2)
  names(weights) <- colnames(R)
  # equity <- matrix(apply(RC, 1, function(x) sum(x * weights)),
  #                  ncol = 1)
  # rownames(equity) <- rownames(R)
  # uvals <- opt$solution[(N + 1):(N + J)]
  # dd <- as.xts(uvals - equity)
  # colnames(dd) <- "DrawDowns"
  # obj <- new("PortMdd", weights = weights, opt = opt, type = "maximum draw-down",
  #            call = call, max_dd = max(dd), DrawDown = dd)
  # return(obj)
  return(weights)
}