R solve.QP not optimal solution mean-variance problem (tracking error minimization)

Question

I have a fictive benchmark of 24 assets named A to X. I have a fictive portfolio which can only trade in a subset (8) of these assets. The lower and upper boundaries indicate my min and max weight in these assets. I want to minimize tracking error of my portfolio to the benchmark. I use solve.QP to get a solution but my solution is not the optimal one.

my cov matrix and boundaries are:

 cov = data.frame(c(0.016,0.013,0.012,0,0.01,0.012,0.003,-0.001,0.002,0.002,0.002,0.002,0.001,0.002,0.001,0.001,0.003,0.002,0,0.001,-0.005,0,0.001,0.002),c(0.013,0.012,0.011,0.001,0.009,0.01,0.003,-0.001,0.002,0.001,0.002,0.002,0.001,0.002,0.001,0,0.002,0.002,0,0.001,-0.004,0,0.001,0.001),c(0.012,0.011,0.011,0.001,0.007,0.011,0.003,0,0.002,0.001,0.002,0.002,0.002,0.001,0.001,0,0.001,0.002,0,0.001,-0.005,0,0.001,0.001),c(0,0.001,0.001,0.012,0.001,0,0.001,-0.001,0.002,-0.001,0.001,0.001,0.001,0,0.001,0.001,-0.002,0.001,0,-0.001,-0.003,-0.001,0.001,-0.002),c(0.01,0.009,0.007,0.001,0.014,0.007,0.004,0,0.002,0.001,0.002,0.003,0,0.005,0.003,0.003,0.003,0.003,-0.002,0.001,-0.005,-0.001,0.001,0.001),c(0.012,0.01,0.011,0,0.007,0.01,0.003,0,0.002,0.001,0.002,0.002,0.002,0.001,0.001,0,0.001,0.002,0.001,0.001,-0.004,0,0.001,0.001),c(0.003,0.003,0.003,0.001,0.004,0.003,0.008,0.001,0.002,0.003,0.002,0.003,0.002,0.001,0.001,0.001,0,0.001,0,-0.001,0,-0.001,0.001,-0.001),c(-0.001,-0.001,0,-0.001,0,0,0.001,0.001,0,0.002,0,0,0,0,0,0,0,0,0,0,0.002,-0.001,0,0),c(0.002,0.002,0.002,0.002,0.002,0.002,0.002,0,0.004,0.001,0.002,0.002,0.002,0.002,0.001,0.001,0,0.001,0.001,0,-0.002,0,0.002,-0.002),c(0.002,0.001,0.001,-0.001,0.001,0.001,0.003,0.002,0.001,0.004,0.001,0.001,0.001,0.001,0,0,0.001,0,0,-0.001,-0.001,-0.001,0,0),c(0.002,0.002,0.002,0.001,0.002,0.002,0.002,0,0.002,0.001,0.002,0.002,0.001,0.002,0.001,0.001,0.001,0.001,0,0,-0.001,0,0,-0.001),c(0.002,0.002,0.002,0.001,0.003,0.002,0.003,0,0.002,0.001,0.002,0.003,0.001,0.002,0.001,0.001,0,0.001,0,0,0,-0.001,0,-0.001),c(0.001,0.001,0.002,0.001,0,0.002,0.002,0,0.002,0.001,0.001,0.001,0.003,0,0,0,0,0,0.001,0,0,0,0.001,-0.001),c(0.002,0.002,0.001,0,0.005,0.001,0.001,0,0.002,0.001,0.002,0.002,0,0.01,0.002,0.002,0.002,0.002,-0.001,0,-0.001,-0.001,0,0),c(0.001,0.001,0.001,0.001,0.003,0.001,0.001,0,0.001,0,0.001,0.001,0,0.002,0.003,0.002,0,0.001,-0.001,0,-0.001,0,0,0),c(0.001,0,0,0.001,0.003,0,0.001,0,0.001,0,0.001,0.001,0,0.002,0.002,0.004,0,0.001,-0.001,0.001,-0.004,0,0,0),c(0.003,0.002,0.001,-0.002,0.003,0.001,0,0,0,0.001,0.001,0,0,0.002,0,0,0.004,0.001,0,0.001,-0.003,0,0,0.001),c(0.002,0.002,0.002,0.001,0.003,0.002,0.001,0,0.001,0,0.001,0.001,0,0.002,0.001,0.001,0.001,0.003,0,0,-0.002,0,0.001,0),c(0,0,0,0,-0.002,0.001,0,0,0.001,0,0,0,0.001,-0.001,-0.001,-0.001,0,0,0.007,0,0.002,0.001,0.001,-0.001),c(0.001,0.001,0.001,-0.001,0.001,0.001,-0.001,0,0,-0.001,0,0,0,0,0,0.001,0.001,0,0,0.003,0.001,0.002,0.002,0),c(-0.005,-0.004,-0.005,-0.003,-0.005,-0.004,0,0.002,-0.002,-0.001,-0.001,0,0,-0.001,-0.001,-0.004,-0.003,-0.002,0.002,0.001,0.042,0.002,0.001,-0.002),c(0,0,0,-0.001,-0.001,0,-0.001,-0.001,0,-0.001,0,-0.001,0,-0.001,0,0,0,0,0.001,0.002,0.002,0.005,0.002,0),c(0.001,0.001,0.001,0.001,0.001,0.001,0.001,0,0.002,0,0,0,0.001,0,0,0,0,0.001,0.001,0.002,0.001,0.002,0.002,-0.001),c(0.002,0.001,0.001,-0.002,0.001,0.001,-0.001,0,-0.002,0,-0.001,-0.001,-0.001,0,0,0,0.001,0,-0.001,0,-0.002,0,-0.001,0.007))

 colnames(cov) = LETTERS[seq( from = 1, to = 24 )]

 rownames(cov) = LETTERS[seq( from = 1, to = 24 )]

 lowerbound = as.matrix(c(0.164,0.131,0.037,0,0.039,0.029,0,0.085,0,0,0.046,0.059,0,0,0,0,0,0,0,0,0,0,0,0))
 upperbound = as.matrix(c(1,1,1,0,1,1,0,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0))

I then optimze as follows:

  require(quadprog)
  library(Matrix)
  Dmat <- nearPD(as.matrix(cov))$mat # make matrix positive definite 
  N = length(cov)
  dvec = rep(0, N)
  A.mat = cbind(rep(1, N), diag(1, N), diag(-1, N))
  bvec = c(1, lowerbound, -upperbound) 

  solution = data.frame(round(solve.QP(as.matrix(Dmat), dvec, A.mat, bvec, meq = 0)$solution,2))
  rownames(solution) = LETTERS[seq( from = 1, to = 24 )]
  solution

This solution is not the optimal solution, asset H gets 50% which is not optimal. Do you know what I am doing wrong?

Update 20200612: I tried a different way using solnp from the Rsolnp package:

  portfolio = c("A", "B", "C", "E", "F", "H", "J", "K") #the assets owned in the portfolio

  weights_bm = c(0.164,0.131,0.037,0.035,0.039,0.029,0.012,0.085,0.014,0.007,0.046,0.059,0.009,0.03,0.019,0.051,0.02,0.015,0.011,0.065, 0.021,0.049,0.047,0.006)

  sol_test <- function(portfolio,weights_bm ,lx,ux,cov){

  cov = as.matrix(rbind(cbind(cov,cov),cbind(cov,cov)))
  wt = rbind(as.matrix(weights_bm*(ux)),as.matrix(-weights_bm))

  obj <- function(x){na.omit(t(x) %*% cov %*% x)} #tracking error

  lx = rbind(as.matrix(lx),as.matrix(wt[(nrow(wt)/2+1):nrow(wt)]))
  ux = rbind(as.matrix(ux),as.matrix(wt[(nrow(wt)/2+1):nrow(wt)]))

  equal <- function(x) {
    equal <- sum(x)
  }

  # solver function
  sol <- solnp(wt,
           fun = obj,
           eqfun = equal,
           eqB = 1,
           LB = as.matrix(lx),
           UB = as.matrix(ux))

   # output from optimalisation
   wt <- sol$pars[1:length(weights_bm )] # optimal portfolio allocation
   names(wt) <- colnames(weights_bm)

   return(wt)

   }

   sol_test(portfolio,weights_bm,lowerbound,upperbound,cov)

However this does not work as the sum of the portfolio weights add up to 0.59 and are just the lowerbound, not optimized outcome.

Answer 1

solved using optimiser from nloptr package