R/weigh_mvef.R

In rengelke/tradr:

#' Title
#'
#' @param R xts object of asset returns
#' @param short argument if short-selling is allowed, "no" by default
#' @param risk_premium risk premium used for modelling
#' @param max_allocation maximum fraction allowed for any of the asset positions
#' @param tangency logical. Should a tangency portfolio, maximum sharpe ratio, be calculated
#'
#' @return xts object containing asset weights
#' @export
#'
#' @import xts
#' @importFrom magrittr %>% %<>%
#'
#' @examples weigh_mvef(stock_returns, max_allocation = 0.05, risk_premium = 0.4)
weigh_mvef <- function (R, short = "no",
                        risk_premium = 0.5,
                        max_allocation = NULL,
                        tangency = FALSE) {

    covariance <- corpcor::make.positive.definite(cov(R))
    n <- ncol(covariance)

    # Create initial amat and bvec assuming only equality constraint
    # (short-selling is allowed, no allocation constraints)
    amat <- matrix (1, nrow = n)
    bvec <- 1
    meq <- 1

    # modify the amat and bvec if short-selling is prohibited
    if(short == "no") {
        amat <- cbind(1, diag(n))
        bvec <- c(bvec, rep(0, n))
    }

    # modify amat and bvec if a max allocation (concentration) is specified
    if(!is.null(max_allocation)) {
        if(max_allocation > 1 | max_allocation < 0){
            stop("max_allocation must be greater than 0 and less than 1")
        }
        if(max_allocation * n < 1) {
            stop("Need to set max_allocation higher; not enough assets to add to 1")
        }
        amat <- cbind(amat, -diag(n))
        bvec <- c(bvec, rep(-max_allocation, n))
    }


    if (tangency) {
        risk_premium_seq <- seq(from = 0, to = 0.5, by = 0.005)
    } else {
        risk_premium_seq <- risk_premium
    }

    eff <- risk_premium_seq %>%
        lapply(function(x) {
            dvec <- colMeans(R) * x
            sol <- quadprog::solve.QP(covariance, dvec = dvec, Amat = amat,
                                      bvec = bvec, meq = meq)
            st_dev <- sqrt(sum(sol$solution * colSums((covariance * sol$solution))))
            exp_ret <- as.numeric(sol$solution %*% colMeans(R))
            sharpe <- exp_ret/st_dev
            c(sol$solution, st_dev, exp_ret, sharpe)
        })

    eff %<>% do.call(rbind, .) %>% as.data.frame() %>%
        `colnames<-`(c(colnames(R), "std_dev", "exp_return", "sharpe"))


    if (tangency) {
        weights <- eff[eff$sharpe == max(eff$sharpe),  ]
    } else {
        weights <- eff
    }

    weights %<>% replace(.<0.0001, 0) %>%
        .[, colnames(.) %in% colnames(R)]
    weights %<>% xts(., order.by = last(index(R)), dateFormat = "Date")

    return(weights)
}