Boston University Questrom School of Business
FE 825 - Spring 2018
Numerical Mean-Variance Optimization in R
Huangyu Chen & Eric Jacquier
The three classic problems in mean-variance portfolio optimization, minimum variance,
maximum expected utility, maximum Sharpe Ratio, can all be written as quadratic pro-
gramming problems. With only expected return, and full investment constraints, they
have analytical solutions. If we add no-short-sales and other bound constraints, one has to
resort to numerical optimization. It is easily done in the EXCEL solver when it is only a
few optimizations. Larger scale projects require programmable languages like R or Matlab.
This note describes the use of solve.QP (package quadprog) and portfolio.optim (pack-
age tseries). It does not repeat the help file you can get from R.
solve.QP(Dmat, dvec, Amat, bvec, meq=0, factorized=FALSE) solves:
min
b
−d′b+ 12b′Db (1)
S.T. A′b ≥ b0 (2)
For portfolio optimization, b is the N-vector of weights, it is an output of the routine, not
an input. D is proportional to the returns covariance matrix. The inputs needed are
• dvec: N-vector, set to zero if not needed.
• Dmat: holds D, a N ×N matrix.
• Amat: holds A, a N ×m matrix of coefficients needed for the m constraints. Note
that solve.QP use A′, not A. Don’t ask why!
• bvec: holds b0 m-vector needed for the m constraints
• meq: number of strict equality constraints. They must be the first (meq) constraints
listed. The remaining (m - meq) constraints are ≥ 0 constraints.
1
1 Minimum variance portfolio
min
w
w′V w
S.T. µ′w = µp
i′w = 1, (3)
where i is a vector of 1s, µp is the desired portfolio mean return (a scalar), and V is the
covariance matrix.
One needs to lump all the constraints in the A′b ≥ b0 formulation needed by solve.QP.
We also input the number of strict equality constraints. We have:
• meq = 2, the mean and full investment constraints are equalities.
• dvec = 0, not needed here
• Dmat = V (or 2V, it does not matter here)
• b0 =
(
µp
1
)
• A′ =
(
µ′
i′
)
, so we have A′w =
(
µp
1
)
.
1.1 No-short-sales constraint
Input them as N additional positivity constraints. Together with the full investment con-
straint, it guarantees that no single weight is above 100% as well.
• Extend the A′ matrix, stacking an identity matrix below the first two rows.
• Extend the b0 vector, adding N zeros.
Write it to see that it does the trick.
1.2 Global MVP
Remove the first constraint (first row of A′ and first element of b0). Don’t forget to change
meq.
1.3 Other constraints
• Incorporate any additional constraint by an additional row to A′ and an additional
element to b0.
• For ≤ constraints, transform them into ≥ constraints by multiplying by −1.
• Put all the equality constraints first and set meq = correctly.
2
1.4 Command portfolio.optim
For the MVP problem subject to a desired return, the command portfolio.optim (package
tseries) is a convenient substitute. It is however a bit limited in scope.
portfolio.optim(x, pm = mean(x), riskless = FALSE, shorts = FALSE, rf =
0.0, reslow = NULL, reshigh = NULL, covmat = cov(x))
• It can’t have sophisticated constraints, but you can directly input bounds on weights
by using shorts = F. You can also use shorts = T. You can alternatively specify
N-vectors of lower bounds and upper bounds for the parameters by setting reslow
and reshigh.
• It computes mean and covariance matrix from the data matrix (input x in the com-
mand). So you can’t use your own vector of means different from the sample mean.
You can however supply an alternative covariance matrix instead of cov(x) by using
covmat = .
• It can’t do a simple global MVP. It always minimizes variance subject to a desired
mean return, pm in the inputs. If you don’t supply pm, it uses the global mean of
all the returns.
• You can set a risk-free rate, using riskless=T, rf = .
2 Maximum Certainty Equivalent Problem
To maximize the CE, with no risk free rate:
max
w
µ′w − 12Aw′V w
S.T. 1′w = 1
Compare to
min
b
−d′b+ 12 b′Db
S.T. A′b ≥ b0
• meq = 1, only the full investment constraints is an equality.
• dvec = µ, the vector of means
• Dmat = γV , where γ is the risk aversion used.
• b0 = 1
• A′ = i′.
• Add the no-short-sales constraint if required.
3
3 Maximum Sharpe Ratio
3.1 Maximum Sharpe Ratio Problem
With the existence of a risk-free asset, the frontier is the CAL. All portfolios on the CAL
have the same Sharpe Ratio as the Tangency Portfolio. Maximizing the Sharpe Ratio
amounts to identifying any one portfolio on the frontier. See the notes. We solve:
min
w
w′V w
S.T (µ−Rf i)′w = µp −Rf ,
where Rf is the risk-free rate.
Choose µp > Rf and use solve.QP with only one equality constraint. Rescale the re-
sulting risky weights so they sum to one to get the Tangency portfolio.
One can also formulate a maximum CE problem with risk-free rate (choose an arbitrary
risk aversion, 2 for example) to find a portfolio on the CAL. The differences with the
previous section are that (1) there is no full investment constraint, and (2) dvec is the
expected return µ−Rf i
4 Beyond
If you’re interested in general numerical optimization, the book Convex Optimization by
Stephen Boyd and Lieven Vandenberghe is a good starting point.
If you need more complex constraints (non-linear constraints) and alternate objective
functions, you will need to use more expansive packages. The two packages below cover a
wide range of portfolio optimization objectives and constraints.
1. Portfolio Analytics
2. fPortfolio
Also, for more optimization package choices in R, take a look at R’s official CRAN Task
page of Optimization and Mathematical Programming.
4
Appendix: solve.QP, Example
\noindent Draw the frontier with 8 industry data .
\begin{lstlisting}[language=R]
require(quadprog)
# Data Preparation and Input Computation
factor <- read.table("factors.txt", header = T)/100
indus<- read.csv("8industries.csv")/100
nasset<- 8
indus <- indus[1:108,] #select time period
Muvec <-apply(indus[,2:9],2,mean)*12
Sigvec<-apply(indus[,2:9],2,sd)*sqrt(12)
Sigmat<-cov(indus[,2:9])*12
ivec<-rep(1,nasset)
\end{lstlisting}
\vspace{0.1in}
\noindent Create vector of $\mu_P$ values, find the corresponding
frontier $\sigma_P$ by solve.QP.
\begin{lstlisting}[language=R]
mufront <- seq(min(Muvec) - 0.1,max(Muvec)+0.15,length = 50)
sdfront <- rep(0,length(mufront))
# formulate standard QP
Dmat <- 2 * Sigmat
dvec <- rep(0, nasset)
Amat <- cbind(Muvec,ivec)
count <- 0
for(mu in mufront){
count <- count + 1
bvec <- rbind(mu,1)
soln <- solve.QP(Dmat, dvec, Amat, bvec, meq=2, factorized=FALSE)
weight <- soln$solution
sdfront[count] <- sqrt(t(weight) %*% Sigmat %*% weight)
}
plot(sdfront,mufront,type = ’l’, col = ’blue’)
\end{lstlisting}
\begin{figure}[htbp]
\begin{center}
\includegraphics[scale = 0.8]{plot1}
\end{center}
5
\end{figure}
\vspace{0.2in}
\begin{lstlisting}[language=R]
# Compare with analytical results
vinv<-solve(Sigmat)
AA <-sum(vinv)
BB <-t(Muvec)%*%vinv%*%ivec
CC <-t(Muvec)%*%vinv%*%Muvec
DEL<-AA*CC-BB^2
sdfront2 <- sqrt((AA*mufront^2-2*BB*mufront + CC)/DEL)
plot(sdfront2,mufront,type = ’l’, col = ’blue’)
\end{lstlisting}
\begin{figure}[htbp]
\begin{center}
\includegraphics[scale = 0.8]{plot2}
\end{center}
\end{figure}
\newpage
\noindent
You could also plot both results on the same plot, to see that they are identical.
\begin{lstlisting}[language=R]
plot(Sigvec,Muvec,type = "p",xlim=c(0.1,0.30),ylim = c(0,0.28),
xlab = "Std. Dev.",ylab = "Mean")
lines(sdfront2,mufront, col = ’red’,lwd=2)
points(sdfront,mufront, col = ’blue’,type=’*’ )
legend(’topleft’,inset = 0.03, legend = c(’quadprog’,’formula’),
lty = c(2,1), col = c(’blue’,’red’))
6








































































































































































































学霸联盟