Why doesn't solve.QP and portfolio.optim gener

Question

The documentation for portfolio.optim {tseries} says that solve.QP {quadprog} is used to generate the solution for finding the tangency portfolio that maximizes the Sharpe ratio. That implies that results should be identical with either function. I'm probably overlooking something, but in this simple example I get similar but not identical solutions for estimating optimal portfolio weights with portfolio.optim and solve.QP. Shouldn't the results be identical? If so, where am I going wrong? Here's the code:

library(tseries)
library(quadprog)


# 1. Generate solution with solve.QP via: comisef.wikidot.com/tutorial:tangencyportfolio

# create artifical data
set.seed(1)
nO     <- 100     # number of observations
nA     <- 10      # number of assets
mData  <- array(rnorm(nO * nA, mean = 0.001, sd = 0.01), dim = c(nO, nA))
rf     <- 0.0001     # riskfree rate (2.5% pa)
mu     <- apply(mData, 2, mean)    # means
mu2    <- mu - rf                  # excess means

# qp
aMat  <- as.matrix(mu2)
bVec  <- 1
zeros <- array(0, dim = c(nA,1))
solQP <- solve.QP(cov(mData), zeros, aMat, bVec, meq = 1)

# rescale variables to obtain weights
w <- as.matrix(solQP$solution/sum(solQP$solution))

# 2. Generate solution with portfolio.optim (using artificial data from above)
port.1 <-portfolio.optim(mData,riskless=rf)
port.1.w <-port.1$pw
port.1.w <-matrix(port.1.w)

# 3. Compare portfolio weights from the two methodologies:

compare <-cbind(w,port.1$pw)

compare

             [,1]        [,2]
 [1,]  0.337932967 0.181547633
 [2,]  0.073836572 0.055100415
 [3,]  0.160612441 0.095800361
 [4,]  0.164491490 0.102811562
 [5,]  0.005034532 0.003214622
 [6,]  0.147473396 0.088792283
 [7,] -0.122882875 0.000000000
 [8,]  0.127924865 0.067705050
 [9,]  0.026626940 0.012507530
[10,]  0.078949672 0.054834759

Answer 1

The difference in your example occurs due to the default value 'shorts=FALSE' in tseries::portfolio.optim(). Therefore you would have to either change the argument or add a non-negativity restriction in your solve.QP problem to reach the same results.

EDIT: While the answer still holds true, there seem to be some other weird default values with tseries::portfolio.optim(). For example it sets the minimum return requirement to pm = mean(x), leading to a random portfolio on the efficiency frontier instead of returning the global minimum variance portfolio if there is no return requirement. Bottom line: Stay with your quadprog::solve.QP solution. Enclosed an example of the wrapper function I use (I just started working with R and while I am quite sure that this delivers mathematically correct results, it might not be the cleanest piece of code):

# --------------------------------------------------------------------------
#' Quadratic Optimization
#' @description Wrapper for quadratic optimization to calculate the general
#'              mean-variance portfolio.
#' @param S           [matrix] Covariance matrix.
#' @param mu          [numeric] Optional. Vector of expected returns.
#' @param wmin        [numeric] Optional. Min weight per asset.
#' @param wmax        [numeric] Optional. Max weight per asset.
#' @param mu_target   [numeric] Optional. Required return, if empty the optimization returns the global minimum variance portfolio
#' @return Returns the mean-variance portfolio or the global minimum variance portfolio
# --------------------------------------------------------------------------

meanvar.pf <- function(S,
                      mu=NULL,
                      wmin=-1000,
                      wmax=1000,
                      mu_target=NULL){
  if (!try(require(quadprog)))
    stop("Execute 'install.packages('quadprog')' and try again")
  if (missing(S))
    stop("Covariance matrix is missing")
  if (!is.null(mu) & dim(S)[1] != length(mu))
    stop("S and mu have non-conformable dimensions")
  N <- ncol(S)
  if (wmin >= 1/N)
    stop("wmin >= 1/N is not feasible")
  if (wmax <= 1/N)
    stop("wmax <= 1/N is not feasible")    
  meq <- 1
  bvec <- c(1, rep(wmin,N), -rep(wmax,N))
  Amat <- cbind(rep(1, N), diag(N), -diag(N))
  if (!is.null(mu_target)) {
    if (is.null(mu))
      stop("Vector of asset returns is missing")
    Amat <- cbind(mu, Amat)
    bvec <- c(mu_target, bvec)
    meq <- 2
  }
  result <- quadprog::solve.QP(Dmat=S, 
                              dvec=rep(0, N), 
                              Amat=Amat, 
                              bvec=bvec, 
                              meq=meq)
  return(result)
}

Answer 2

The one and the only way to deal with such situations is to browse the source. In your case, it is accessible via tseries:::portfolio.optim.default.

Now, to find the difference between those two calls, we may narrow down the issue by defining an equivalent helper function:

foo <- function(x, pm = mean(x), covmat = cov(x), riskless = FALSE, rf = 0) 
{
  x <- mData
  pm <- mean(x)
  covmat <- cov(x)
  k <- dim(x)[2]
  Dmat <- covmat
  dvec <- rep.int(0, k)
  a1 <- colMeans(x) - rf
  a2 <- matrix(0, k, k)
  diag(a2) <- 1
  b2 <- rep.int(0, k)
  Amat <- t(rbind(a1, a2))
  b0 <- c(pm - rf, b2)
  solve.QP(Dmat, dvec, Amat, bvec = b0, meq = 1)$sol
}

identical(portfolio.optim(mData, riskless=TRUE, rf=rf)$pw,
          foo(mData, riskless=TRUE, rf=rf))
#[1] TRUE

With that, you can see that 1) riskless=rf is not the intended way, riskless=TRUE, rf=rf is the correct one; 2) there are several differences in Amat and bvec.

I am not an expert in portfolio optimization, so I do not know what's the explanation behind these additional constraints and if they should be there in the first place, but at least you can see what exactly causes the difference.