---
title: "Portfolio Theory with No-Short Sales"
author: "Eric Zivot"
date: "August 11, 2015"
output: slidy_presentation
---

## Set options and load packages

```{r}
options(digits=3, width=70)
# install IntroCompFinR package from R-forge
# use install.packages("IntroCompFinR", repos="http://R-Forge.R-project.org")
library(quadprog)
library(IntroCompFinR)
library(PerformanceAnalytics)
library(zoo)
Sys.setenv(TZ="UTC")
cex.val = 2
```

## Three asset example data

Estimates of CER model for Microsoft, Nordstrom and Starbucks stock from monthly returns over the period January 1995 to January 2000.

```{r}
asset.names <- c("MSFT", "NORD", "SBUX")
mu.vec = c(0.0427, 0.0015, 0.0285)
names(mu.vec) = asset.names
sigma.mat = matrix(c(0.0100, 0.0018, 0.0011,
		   0.0018, 0.0109, 0.0026,
		   0.0011, 0.0026, 0.0199),
		 nrow=3, ncol=3)
dimnames(sigma.mat) = list(asset.names, asset.names)
r.f = 0.005
mu.vec
sd.vec = sqrt(diag(sigma.mat))
cov2cor(sigma.mat)
```

## Risk return characteristics

```{r, echo=FALSE}
plot(sd.vec, mu.vec,  ylim=c(0, 0.06), xlim=c(0, 0.20), ylab=expression(mu[p]),
     xlab=expression(sigma[p]), pch=16, col="blue", cex=2.5, cex.lab=1.75)     
text(sd.vec, mu.vec, labels=asset.names, pos=4, cex = cex.val)
text(0, r.f, labels=expression(r[f]), pos=2, cex = cex.val)

```

## Random portfolios: no short sales

```{r}
set.seed(123)
x.msft = runif(400, min=0, max=1)
x.nord = runif(400, min=0, max=1)
x.sbux = 1 - x.msft - x.nord
long.only = which(x.msft > 0 & x.nord > 0 & x.sbux > 0)
x.msft = x.msft[long.only]
x.nord = x.nord[long.only]
x.sbux = x.sbux[long.only]
length(long.only)
```

## Random portfolios: no short sales

```{r}
# plot portfolio weights
chart.StackedBar(cbind(x.msft, x.nord, x.sbux), 
                 main="Random long-only portfolio weight vectors",
                 xlab = "portfolio", ylab = "weights",
                 xaxis.labels=as.character(1:length(long.only)))
```

## Risk-return characteristics: long-only random portfolios

```{r, echo=FALSE}
plot(sd.vec, mu.vec,  ylim=c(0, 0.05), xlim=c(0, 0.2), ylab=expression(mu[p]),
     xlab=expression(sigma[p]), type="n")     
for (i in 1:length(long.only)) {
  z.vec = c(x.msft[i], x.nord[i], x.sbux[i])
  mu.p = crossprod(z.vec,mu.vec)
  sig.p = sqrt(t(z.vec)%*%sigma.mat%*%z.vec)
  points(sig.p, mu.p, pch=16, col="grey", cex=1.5)
}
points(sd.vec, mu.vec, pch=16, col="blue", cex=2.5, cex.lab=1.75)
text(sd.vec, mu.vec, labels=asset.names, pos=4, cex = cex.val)
```


## Global minimum variance portfolio: no constraints on short sales

Use the **IntroCompFinR** package function <code>globalMin.portfolio()</code>

```{r}
gmin.port = globalMin.portfolio(mu.vec, sigma.mat)
gmin.port

```

* Global minimum variance portfolio without contraints on short sales does not have any negative weights (shorts)
* Short sales constraint is not binding

## Use the quadprog function <code>solve.QP()</code> to compute global minimum variance portfolio

First we set up the appropriate restriction matrices required by the function <code>slove.QP()</code>

```{r}
D.mat = 2*sigma.mat
d.vec = rep(0, 3)
A.mat = cbind(rep(1,3), diag(3))
b.vec = c(1, rep(0,3))
```

```{r}
D.mat
```

```{r}
d.vec
```

```{r}
t(A.mat)
```

```{r}
b.vec
```

## Use the quadprog function <code>solve.QP()</code> to compute global minimum variance portfolio

The function <code>solve.QP()</code> takes the restrictions matrices as inputs

```{r}
args(solve.QP)
```

<code>solve.QP()</code> returns a list containing information about the optimization

```{r}
qp.out = solve.QP(Dmat=D.mat, dvec=d.vec,
                  Amat=A.mat, bvec=b.vec, meq=1)
class(qp.out)
```

The components of the returned list are

```{r}
names(qp.out)
```

The portfolio weights are in the <code>solution</code> component

```{r}
qp.out$solution
```

The solution satisfies the required constraints (weights sum to one and all weights are positive)

```{r}
sum(qp.out$solution)
```

The minimized value of the objective function (portfolio variance) is in the <code>value</code> component

```{r}
qp.out$value
```

## Use the quadprog function <code>solve.QP()</code> to compute global minimum variance portfolio

Compute the mean and volatility of the minimum variance portfolio

```{r}
# Expected return
er.gmin.ns = as.numeric(crossprod(qp.out$solution, mu.vec))
# Volatility
sd.gmin.ns = sqrt(qp.out$value)
```

```{r}
er.gmin.ns
```

```{r}
sd.gmin.ns
```

## Use the IntroCompFinR function <code>globalMin.Portfolio()</code> to compute the minimum variance portfolio

Use the optional argument <code>shorts=TRUE</code> to impose no-short sales restrictions

```{r}
args(globalMin.portfolio)
```

The function sets up the restriction matrices for <code>solve.QP()</code> to compute the global minimum variance portfolio.

```{r}
gmin.port = globalMin.portfolio(mu.vec, sigma.mat,
                                shorts=FALSE)
gmin.port
```

## Compute minimum variance portfolio with same mean as Microsoft

First find the minimum variance portfolio allowing short sales

```{r}
eMsft.port = efficient.portfolio(mu.vec, sigma.mat, target.return = mu.vec["MSFT"])
eMsft.port
```


## Compute minimum variance portfolio with same mean as Microsoft

First we set up the appropriate restriction matrices required by the function <code>slove.QP()</code>

```{r}
D.mat = 2*sigma.mat
d.vec = rep(0, 3)
A.mat = cbind(mu.vec, rep(1,3), diag(3))
b.vec = c(mu.vec["MSFT"], 1, rep(0,3))
```


```{r}
D.mat
```

```{r}
d.vec
```

```{r}
A.mat
```

```{r}
b.vec
```

## Compute minimum variance portfolio with same mean as Microsoft

Next we use <code>solve.QP()</code> to find the solution

```{r}
qp.out = solve.QP(Dmat=D.mat, dvec=d.vec,
                  Amat=A.mat, bvec=b.vec, meq=2)
names(qp.out$solution) = names(mu.vec)
round(qp.out$solution, digits=3)
```

* The argument <code>meq = 2</code> specifies 2 equality constraints. It tells <code>solve.QP()</code> how to determine $A_{eq}$ and $A_{neq}$ from $A$

* Notice that the no-short sales solution forces the weights on Nordstrom and Starbucks to 0

* It is not possible to find an no-short sale efficient portfolio that has a higher mean than Microsoft

## Compute minimum variance portfolio with same mean as Microsoft

The IntroCompFinR function <code>efficient.portfolio()</code> with <code>shorts=FALSE</code> uses <code>solve.QP()</code> to compute a no-short sales efficient portfolio

```{r}
efficient.portfolio(mu.vec, sigma.mat, target.return = mu.vec["MSFT"], shorts = FALSE)
```


## Infeasible no-short sale efficient portfolio

Suppose you try to find an efficient portfolio with target return higher than the mean for Microsoft

```{r}
# efficient.portfolio(mu.vec, sigma.mat, target.return = mu.vec["MSFT"]+0.01, shorts = FALSE)
# Error in quadprog::solve.QP(Dmat = Dmat, dvec = dvec, Amat = Amat, bvec = bvec,  : 
#  constraints are inconsistent, no solution!
```

* Notice that <code>solve.QP()</code> throws an error saying there is no compatible solution

## Compute efficient frontier allowing for short sales

Use the IntroCompFinR function <code>efficient.frontier()</code> to compute the efficient set of frontier portfolios that allow short sales

```{r}
ef = efficient.frontier(mu.vec, sigma.mat, alpha.min=0, 
                         alpha.max=1, nport=10)
ef$weights

```

* Notice that two efficient portfolios have negative weights (shorts)
* Portfolio 10 is the global minimum variance portfolio

## Compute efficient frontier allowing for short sales

```{r, echo=FALSE}
plot(ef$sd, ef$er, type="b", ylim=c(0.02, 0.05), xlim=c(0.06, 0.11), 
     pch=16, col="blue", cex=2, ylab=expression(mu[p]), xlab=expression(sigma[p]))
points(sd.vec, mu.vec, pch=16, col="black", cex=2)
text(sd.vec, mu.vec, labels=asset.names, pos=4)
text(ef$sd[1], ef$er[1], labels="port 10", pos=4)
text(ef$sd[2], ef$er[2], labels="port 9", pos=4)
```

* port 1 and port 2 have negative weights in Nordstrom

## Compute efficient frontier without short sales

It is easy to compute the efficient frontier without short sales by running a simple loop

```{r}
mu.vals = seq(gmin.port$er, max(mu.vec), length.out=10)
w.mat = matrix(0, length(mu.vals), 3)
sd.vals = rep(0, length(sd.vec))
colnames(w.mat) = names(mu.vec)
D.mat = 2*sigma.mat
d.vec = rep(0, 3)
A.mat = cbind(mu.vec, rep(1,3), diag(3))
for (i in 1:length(mu.vals)) {
  b.vec = c(mu.vals[i],1,rep(0,3))
  qp.out = solve.QP(Dmat=D.mat, dvec=d.vec,
                    Amat=A.mat, bvec=b.vec, meq=2)
  w.mat[i, ] = qp.out$solution
  sd.vals[i] = sqrt(qp.out$value)
}
```

## Compare frontiers

```{r, echo=FALSE}
plot(ef$sd, ef$er, type="b", ylim=c(0.02, 0.05), xlim=c(0.06, 0.11), 
     pch=16, col="blue", cex=2, ylab=expression(mu[p]), xlab=expression(sigma[p]))
points(sd.vec, mu.vec, pch=16, col="black", cex=2)
text(sd.vec, mu.vec, labels=asset.names, pos=4)
points(sd.vals, mu.vals, type="b", pch=16, col="red", cex=1.5)
text(ef$sd[1], ef$er[1], labels="port 10", pos=2)
text(ef$sd[2], ef$er[2], labels="port 9", pos=2)
```

* Here, the no-short sales efficient frontier lies "inside" and "to the right" of the short sales frontier

## Compute efficient frontier without short sales

```{r}
ef.ns = efficient.frontier(mu.vec, sigma.mat, alpha.min=0, 
                         alpha.max=1, nport=10, shorts=FALSE)
ef.ns$weights
```

## Copmpare efficient frontiers

```{r}
plot(ef$sd, ef$er, type="b", ylim=c(0.001, 0.05), xlim=c(0.06, 0.15), 
     pch=16, col="blue", cex=2, ylab=expression(mu[p]), xlab=expression(sigma[p]))
points(sd.vec, mu.vec, pch=16, col="black", cex=2)
text(sd.vec, mu.vec, labels=asset.names, pos=4)
points(ef.ns$sd, ef.ns$er, type="b", pch=16, col="red")
text(ef$sd[1], ef$er[1], labels="port 10", col="blue", pos=3)
text(ef$sd[2], ef$er[2], labels="port 9", col="blue", pos=3)
text(ef.ns$sd[10], ef.ns$er[10], labels="port 10", col="red", pos=1)
text(ef.ns$sd[9], ef.ns$er[9], labels="port 9", col="red", pos=1)
```


