xuebaunion@vip.163.com

3551 Trousdale Rkwy, University Park, Los Angeles, CA

留学生论文指导和课程辅导

无忧GPA：https://www.essaygpa.com

工作时间：全年无休-早上8点到凌晨3点

微信客服：xiaoxionga100

微信客服：ITCS521

r代写-FE 825

时间：2021-03-13

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

学霸联盟

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

学霸联盟