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
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.