R/optimalPortfolio.R

In ExpectedReturns: Computation of implied expected returns

optimalPortfolio = function(Sigma, mu = NULL, semiDev = NULL, control = list()){ 
  #########################################################################################
  # Computes the weight of the chosen portfolio
  # INPUTs
  #   Sigma        : matrix (N x N) covariance matrix
  #   mu           : vector (N x 1) expected returns (optional)
  #   semiDev      : vector (N x 1) semi-deviation (optional)
  #   control      : a control list 
  #   The argument control is a list that can supply any of the following components
  #     type       : "mv", "minvol", "erc", "maxdiv", "riskeff"
  #     constraint : "none", "lo" (long only), "gross"
  #     gross.c    : default = 1.6 (130%-30%)
  #     gamma      : risk-averesion default = 0.89
  # OUTPUTs
  #   w            : vector (N x 1) optimal weights
  ######################################################################################### 
  if (missing(Sigma)){
    stop ('A covariance matrix (Sigma) is required')
  }
  if (!is.matrix(Sigma)){
    stop ('Sigma must be a matrix')
  }
#   if (!identical(t(Sigma), Sigma)){
#     stop ('Sigma must be a symmetric matrix')
#   }
  if (!isSymmetric(Sigma)){
    stop ('Sigma must be a symmetric matrix')
  }  

  ctr = .ctrPortfolio(control)
  
  if (ctr$type[1] == "mv"){
    w = .mvPortfolio(mu = mu, Sigma = Sigma, control = control)
  }
  else if (ctr$type[1] == "minvol"){
    w = .minvolPortfolio(Sigma = Sigma, control = control)
  }
  else if (ctr$type[1] == "erc"){
    w = .ercPortfolio(Sigma = Sigma, control = control)
  }
  else if (ctr$type[1] == "maxdiv"){
    w = .maxdivPortfolio(Sigma = Sigma, control = control)
  }
  else if (ctr$type[1] == "riskeff"){
    w = .riskeffPortfolio(Sigma = Sigma, semiDev = semiDev, control = control)
  }
  else{
    stop('control$type is not well defined')
  }
  return(w)
}

.ctrPortfolio = function(control = list()){
  #########################################################################################  
  # Function used to control the list input
  # INPUTs
  #   control : a control list 
  #   The argument control is a list that can supply any of the following components
  #     type : "mv", "minvol", "erc", "maxdiv", "riskeff"
  #     constraint : "none", "lo", "gross"
  #     gross.c    : default = 1.6
  #     gamma      : default = 0.89
  # OUTPUTs
  #   control : list
  #########################################################################################  
  if (!is.list(control)){
    stop ('control must be a list') 
  }
  if (length(control) == 0){
    control = list(type = "mv", constraint = "none", gross.c = 1.6, LB = NULL, UB = NULL, gamma = 0.89)
  }
  nam = names(control)
  if (!("type" %in% nam) || is.null(control$type)){
    control$type = c("mv", "minvol", "erc", "maxdiv", "riskeff")
  }
  if (!("constraint" %in% nam) || is.null(control$constraint)){
    control$constraint = c("none", "lo", "gross")
  }
  if (!("gross.c" %in% nam) || is.null(control$gross.c)){
    control$gross.c = 1.6
  }
  if (!("LB" %in% nam) || is.null(control$LB)){
    control$LB = NULL
  }
  if (!("UB" %in% nam) || is.null(control$UB)){
    control$UB = NULL
  }
  if (!("gamma" %in% nam) || is.null(control$gamma)){
    control$gamma = c(0.8773, 2.7063, 3.7950)
  }
  return(control)
}

.mvPortfolio = function(mu, Sigma, control = list()){
  #########################################################################################  
  # Compute the weight of the mean-variance portfolio
  # INPUTs
  #   Sigma : matrix (N x N) covariance matrix
  # OUTPUTs
  #   w     : vector (N x 1) weight
  #########################################################################################    
  ctr = .ctrPortfolio(control)
  
  if (is.null(mu)){
    stop ('A vector of mean (mu) is required to compute the mean-variance portfolio')
  }
  if (ctr$constraint[1] == "none"){
    invSigmamu = solve(Sigma, mu)
    w = (1 / ctr$gamma[1]) * invSigmamu / sum(invSigmamu)  # David, rf = 0 !??
  }
  else if (ctr$constraint[1] == "lo"){
    n = dim(Sigma)[1]
    Dmat = ctr$gamma[1] * Sigma
    Amat = cbind(rep(1,n), diag(n))
    bvec = c(1, rep(0,n))
    w    = solve.QP(Dmat = Dmat, dvec = mu, Amat = Amat, bvec = bvec, meq = 1)$solution   
  }
  else if (ctr$constraint[1] == "gross"){
    
    .meanvar = function(w){
      Sigmaw = crossprod(Sigma, w)
      opt    = - as.numeric( crossprod(mu, w) ) + 0.5 * ctr$gamma[1] * as.numeric( crossprod(w, Sigmaw) ) 
      return(opt)
    }
    
    n  = dim(Sigma)[1]
    w = rep(1/n, n)  
    w = slsqp(x0 = w, fn = .meanvar, hin = .grossConstraint, heq = .eqConstraint, nl.info = FALSE, control = list(xtol_rel = 1e-18, check_derivatives = FALSE))$par
    
    
#     n  = dim(Sigma)[1]
#     
#     Aeq = c(rep(1,n),rep(0,n))
#     beq = 1
#     
#     A = rbind( Aeq, cbind(diag(n), -diag(n)), cbind(-diag(n), -diag(n)), c(rep(0,n),rep(1,n)) )
#     bupper = c(beq, rep(0,2*n), 1.6)
#     blower = c(beq, rep(-Inf, 2*n+1))
#     
#     .meanvar = function(w){
#       w = w[1:n]
#       Sigmaw = crossprod(Sigma, w)
#       opt    = - as.numeric( crossprod(mu, w) ) + 0.5 * ctr$gamma[1] * as.numeric( crossprod(w, Sigmaw) ) 
#       return(opt)
#     }
#     
#     n  = dim(Sigma)[1]
#     par = rep(1/n, 2*n)  
#     nlin = list()
#     
#     w = donlp2( 
#       par = par, fn = .meanvar,
#       par.upper = rep(+Inf, length(par)), 
#       par.lower = rep(-Inf, length(par)),   
#       A = A,
#       lin.upper = bupper, 
#       lin.lower = blower, 
#       nlin = list(),
#       nlin.upper = rep(+Inf, length(nlin)), 
#       nlin.lower = rep(-Inf, length(nlin)),
#       control = donlp2Control(iterma = 200000),
#       control.fun = function(lst){return(TRUE)},
#       env = .GlobalEnv, 
#       name = NULL)$par
#     w = w[1:n]
    
  }
  else{
    stop ('control$constraint not well defined')
  }
  return(w)
}

.minvolPortfolio = function(Sigma, control = list() ){
  #########################################################################################  
  # Compute the weight of the minimum volatility portfolio
  # INPUTs
  #   Sigma : matrix (N x N) covariance matrix
  # OUTPUTs
  #   w     : vector (N x 1) weight
  #########################################################################################   
  ctr = .ctrPortfolio(control)
  n = dim(Sigma)[1]
  
  if (ctr$constraint[1] == "none"){
    
    tmp = solve(Sigma, rep(1, n))
    w   = tmp / sum(tmp)
    #dvec = vector("double", n)
    #Amat = as.matrix(rep(1,n))
    #bvec = 1
    #solve.QP(Dmat = Sigma, dvec = dvec, Amat = Amat, bvec = bvec, meq = 1)
    
  }
  else if (ctr$constraint[1] == "lo"){
    dvec = rep(0, n)
    Amat = cbind(rep(1, n), diag(n))
    bvec = c(1, rep(0, n))
    w = solve.QP(Dmat = Sigma, dvec = dvec, Amat = Amat, bvec = bvec, meq = 1)$solution    
  }
  else if (ctr$constraint[1] == "gross"){
    
    .minvol = function(w){
      Sigmaw = crossprod(Sigma, w)
      v = as.numeric(crossprod(w, Sigmaw))
      return(v)
    }
    
    n  = dim(Sigma)[1]
    w = rep(1/n, n)  
    w = slsqp(x0 = w, fn = .minvol, hin = .grossConstraint, heq = .eqConstraint, nl.info = FALSE, control = list(xtol_rel = 1e-18, check_derivatives = FALSE))$par
     
#     w = rep(1/n, n) 
#     
#     .minv = function(w){
#       w = w[1:n]
#       Sigmaw = crossprod(Sigma, w)
#       v = as.numeric(crossprod(w, Sigmaw))
#       return(v)
#     }
#     
#     Aeq = c(rep(1,n),rep(0,n))
#     beq = 1
#     
#     A = rbind( Aeq, cbind(diag(n), -diag(n)), cbind(-diag(n), -diag(n)), c(rep(0,n),rep(1,n)) )
#     bupper = c(beq, rep(0,2*n), 1.6)
#     blower = c(beq, rep(-Inf, 2*n+1))
#     
#     n  = dim(Sigma)[1]
#     par = rep(1/n, 2*n)  
#     nlin = list()
#     
#     w = donlp2( 
#       par = par, fn = .minv,
#       par.upper = rep(+Inf, length(par)), 
#       par.lower = rep(-Inf, length(par)),   
#       A = A,
#       lin.upper = bupper, 
#       lin.lower = blower, 
#       nlin = list(),
#       nlin.upper = rep(+Inf, length(nlin)), 
#       nlin.lower = rep(-Inf, length(nlin)),
#       control = donlp2Control(iterma = 200000),
#       control.fun = function(lst){return(TRUE)},
#       env = .GlobalEnv, 
#       name = NULL)$par
#     w = w[1:n]

  }
  else{    
    stop ('control$constraint not well defined')    
  }
  return(w)  
}

.ercPortfolio = function(Sigma, control = list()){
  #########################################################################################  
  # Compute the weight of the equal-risk-contribution portfolio
  # INPUTs
  #   Sigma : matrix (N x N) covariance matrix
  # OUTPUTs
  #   w     : vector (N x 1) weight
  #########################################################################################   
  ctr = .ctrPortfolio(control)
  
  n = dim(Sigma)[2]
  w = rep(1/n, n)  
  
  .pRC = function(w){
    Sigmaw = crossprod(Sigma, w)
    pRC    = (w * Sigmaw) / as.numeric( crossprod(w, Sigmaw) )
    d = sum( (pRC - 1/n)^2 )
    return(d)
  }
  
  if (ctr$constraint[1] == "none"){
    w = slsqp(x0 = w, fn = .pRC, heq = .eqConstraint, lower = rep(0,n), nl.info = FALSE, control = list(xtol_rel = 1e-8, check_derivatives = FALSE))$par  
  }
  else if(ctr$constraint[1] == "lo"){ 
    w = slsqp(x0 = w, fn = .pRC, heq = .eqConstraint, lower = rep(0,n), nl.info = FALSE, control = list(xtol_rel = 1e-8, check_derivatives = FALSE, maxeval = 2000))$par  
  }
  else if(ctr$constraint[1] == "gross"){    
    w = slsqp(x0 = w, fn = .pRC, heq = .eqConstraint, lower = rep(0,n),  nl.info = FALSE, control = list(xtol_rel = 1e-8, check_derivatives = FALSE, maxeval = 2000))$par     
  }
  else{
    stop ('control$constraint not well defined')
  }
  return(w)
}

.maxdivPortfolio = function(Sigma, control = list()){ 
  #########################################################################################  
  # Compute the weight of weight of the maximum diversification portfolio
  # INPUTs
  #   Sigma : matrix (N x N) covariance matrix
  # OUTPUTs
  #   w     : vector (N x 1) weight
  ######################################################################################### 
  ctr = .ctrPortfolio(control)
  
  n = dim(Sigma)[2]
  w = rep(1/n, n)
  
  .divRatio = function(w){
    sig      = sqrt(diag(Sigma))
    Sigmaw   = crossprod(Sigma, w)
    divRatio = as.numeric(- crossprod(w, sig) / sqrt( crossprod(w, Sigmaw)) )
    return(divRatio)
  }
  if (ctr$constraint[1] == "none"){
    w = slsqp(x0 = w, fn = .divRatio, heq = .eqConstraint, nl.info = FALSE, control = list(xtol_rel = 1e-8, check_derivatives = FALSE, maxeval = 2000))$par
  }
  else if(ctr$constraint[1] == "lo"){ 
    w = slsqp(x0 = w, fn = .divRatio, lower = rep(0,n), heq = .eqConstraint, nl.info = FALSE, control = list(xtol_rel = 1e-8, check_derivatives = FALSE, maxeval = 2000))$par    
  }
  else if(ctr$constraint[1] == "gross"){
    w = slsqp(x0 = w, fn = .divRatio, hin = .grossConstraint, heq = .eqConstraint, nl.info = FALSE, control = list(xtol_rel = 1e-8, check_derivatives = FALSE, maxeval = 2000))$par
  }
  else{
    stop ('control$constraint not well defined')
  }
  return(w)
}

.riskeffPortfolio = function(Sigma, semiDev, control = list()){
  #########################################################################################  
  # Compute the weight of the risk-efficient portfolio
  # INPUTs
  #   Sigma : matrix (N x N) covariance matrix
  # OUTPUTs
  #   w     : vector (N x 1) weight
  ######################################################################################### 
  ctr = .ctrPortfolio(control)
  
  if (is.null(semiDev)){
    stop ('A vector of semideviation (semiDev) is require to compute the risk-efficient portfolio')
  }
  
  n       = dim(Sigma)[2]
  pct     = c(0, quantile(semiDev, probs = seq(0.1, 1, 0.1))) 
  epsilon = vector("double", n)
  J       = matrix(rep(0, n^2), ncol = n)  
  
  for (i in 2:11) {
    pos = semiDev > pct[i-1] & semiDev <= pct[i]
    J[pos, i - 1] = 1
    epsilon[i - 1] = median(semiDev[pos])
  }
  Jepsilon = crossprod(t(J), epsilon)
  # !!! TOFIX !!! additional constraints used to stabilize optimization
  LB = (1 / (2 * n)) * rep(1, n)
  UB = (2 / n) * rep(1, n)
  
  w = (UB - LB)
  w = w / sum(w)
  
  .distRiskEff = function(w){
    Sigmaw = crossprod(Sigma, w)
    d      = as.numeric(-crossprod(w, Jepsilon) / sqrt( crossprod(w, Sigmaw)))
    return(d)
  }
  
  .eqConstraintReff = function(w){
    return(sum(w))
  }
  
  if (ctr$constraint[1] == "none"){
    w = slsqp(x0 = w, fn = .distRiskEff, heq = .eqConstraint, lower = LB, upper = UB, nl.info = FALSE, control = list(xtol_rel = 1e-8, check_derivatives = FALSE, maxeval = 2000))$par
  }
  else if(ctr$constraint[1] == "lo"){ 
    w = slsqp(x0 = w, fn = .distRiskEff, heq = .eqConstraint, lower = LB, upper = UB, nl.info = FALSE, control = list(xtol_rel = 1e-8, check_derivatives = FALSE, maxeval = 2000))$par    
  }
  else if(ctr$constraint[1] == "gross"){
    w = slsqp(x0 = w, fn = .distRiskEff, heq = .eqConstraint, lower = LB, upper = UB, nl.info = FALSE, control = list(xtol_rel = 1e-8, check_derivatives = FALSE, maxeval = 2000))$par
  }
  else{
    stop ('control$constraint not well defined')
  }
  return(w)
}

#########################################################################################  
## Constraints used by the optimizers
#########################################################################################  
.eqConstraint = function(w){
  return(sum(w) - 1)
}

.grossConstraint = function(w){
  return(1.6 - norm(as.matrix(w), type = "1") )
}