---
title: "Portfolio Theory with Matrix Algebra"
author: "Eric Zivot"
date: "August 6, 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(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
cov2cor(sigma.mat)
```

## Risk return characteristics

```{r, echo=FALSE}
sd.vec = sqrt(diag(sigma.mat))
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)

```

* MSFT has highest Sharpe ratio
* NORD has lowest Sharpe ratio

## Example portfolio: equally weighted

```{r}
# create vector of portfolio weights
x.vec = rep(1,3)/3
names(x.vec) = asset.names
# compute mean, variance and volatility
mu.p.x = crossprod(x.vec,mu.vec)
sig2.p.x = t(x.vec)%*%sigma.mat%*%x.vec
sig.p.x = sqrt(sig2.p.x)
```

```{r}
# show mean and volatility
mu.p.x
sig.p.x

```

## Example portfolio: long-short

```{r}
# create vector of portfolio weights
y.vec = c(0.8, 0.4, -0.2)
names(y.vec) = asset.names
# compute mean, variance and volatility
mu.p.y = crossprod(y.vec,mu.vec)
sig2.p.y = t(y.vec)%*%sigma.mat%*%y.vec
sig.p.y = sqrt(sig2.p.y)
```

```{r}
# show mean and volatility
mu.p.y
sig.p.y
```

## Covariance and correlation between example portfolio returns

```{r}
# covariance
sig.xy = t(x.vec)%*%sigma.mat%*%y.vec
sig.xy
```

```{r}
# correlation
rho.xy = sig.xy/(sig.p.x*sig.p.y)
rho.xy

```

## Risk-return characteristics: example portfolios

```{r, echo=FALSE}
cex.val = 1.15
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)
points(sig.p.x, mu.p.x, pch=16, col="black", cex=2.5)
text(sig.p.x, mu.p.x, labels="EQUAL WEIGHT", pos=4, cex = cex.val)
points(sig.p.y, mu.p.y, pch=16, col="black", cex=2.5)
text(sig.p.y, mu.p.y, labels="LONG-SHORT", pos=4, cex = cex.val)
text(0, r.f, labels=expression(r[f]), pos=2, cex = cex.val)

```

## Risk-return characteristics: random portfolios

Create 100 random portfolio vectors with weights that sum to one.

```{r}
set.seed(123)
x.msft = runif(100, min=-1.5, max=1.5)
x.nord = runif(100, min=-1.5, max=1.5)
x.sbux = 1 - x.msft - x.nord
head(cbind(x.msft, x.nord, x.sbux))
```

## Risk-return characteristics: random portfolios

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

## Risk-return characteristics: random portfolios

```{r, echo=FALSE}
plot(sd.vec, mu.vec,  ylim=c(-0.03, 0.08), xlim=c(0, 0.4), 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)
for (i in 1:length(x.msft)) {
  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="black", cex=1.5)
}

```

* With more than two assets, set of feasible portfolios is no longer one side of a hyperbole
* Set of feasible portfolios is a solid space
* Efficient portfolios are on the upper boundary (above minimimum variance portfolio)

## Compute global minimum variance portfolio: Method 1

Use $z_m = A^{-1}_m b$.

```{r}
top.mat = cbind(2*sigma.mat, rep(1, 3))
bot.vec = c(rep(1, 3), 0)
Am.mat = rbind(top.mat, bot.vec)
b.vec = c(rep(0, 3), 1)
z.m.mat = solve(Am.mat)%*%b.vec
m.vec = z.m.mat[1:3,1]
# minimum variance portfolio weights
m.vec
```

Mean and volatility of minimum variance portfolio

```{r}
mu.gmin = as.numeric(crossprod(m.vec, mu.vec))
sig2.gmin = as.numeric(t(m.vec)%*%sigma.mat%*%m.vec)
sig.gmin = sqrt(sig2.gmin)
mu.gmin
sig.gmin
```

## Compute global minimum variance portfolio: Method 2

Use analytic formula for minimum variance portfolio

```{r}
one.vec = rep(1, 3)
sigma.inv.mat = solve(sigma.mat)
top.mat = sigma.inv.mat%*%one.vec
bot.val = as.numeric((t(one.vec)%*%sigma.inv.mat%*%one.vec))
m.mat = top.mat/bot.val
m.mat[, 1]

```

## Show minimum variance portfolio 

```{r, echo=FALSE}
plot(sd.vec, mu.vec,  ylim=c(-0.03, 0.08), xlim=c(0, 0.4), 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)
points(sig.gmin, mu.gmin, pch=16, cex=2.5, col="green")
text(sig.gmin, mu.gmin, labels="GLOBAL MIN", pos=2.5, cex = cex.val)
text(0, r.f, labels=expression(r[f]), pos=2, cex = cex.val)
for (i in 1:length(x.msft)) {
  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="black", cex=1.5)
}


```

* If there is "justice in the world" then the mean and volatility of the global minimum variance portfolio will plot at the tip of the Markowitz bullet.
* It does!

## Compute efficient portfolio with the same mean as Microsoft

Use matrix algebra formula to compute efficient portfolio.

```{r}
top.mat = cbind(2*sigma.mat, mu.vec, rep(1, 3))
mid.vec = c(mu.vec, 0, 0)
bot.vec = c(rep(1, 3), 0, 0)
Ax.mat = rbind(top.mat, mid.vec, bot.vec)
bmsft.vec = c(rep(0, 3), mu.vec["MSFT"], 1)
z.mat = solve(Ax.mat)%*%bmsft.vec
x.vec = z.mat[1:3,]
x.vec
```

Compute mean and volatility of efficient portfolio.

```{r}
mu.px = as.numeric(crossprod(x.vec, mu.vec))
sig2.px = as.numeric(t(x.vec)%*%sigma.mat%*%x.vec)
sig.px = sqrt(sig2.px)
mu.px
sig.px
```

Compare with mean and volatility of Microsoft.

```{r}
mu.vec["MSFT"]
sd.vec["MSFT"]
```

* Efficient portfolio has slightly smaller volatility than Microsoft.
* Microsoft is near the efficient frontier boundary

## Compute efficient portfolio with the same mean as Starbucks

```{r}
# solve for portfolio weights
bsbux.vec = c(rep(0, 3), mu.vec["SBUX"], 1)
z.mat = solve(Ax.mat)%*%bsbux.vec
y.vec = z.mat[1:3,]
y.vec
```

```{r}
# compute mean, variance and std deviation
mu.py = as.numeric(crossprod(y.vec, mu.vec))
sig2.py = as.numeric(t(y.vec)%*%sigma.mat%*%y.vec)
sig.py = sqrt(sig2.py)
mu.py
sig.py
```

```{r}
# compare with Starbucks
mu.vec["SBUX"]
sd.vec["SBUX"]
```

* Efficient portfolio has much smaller volatility than Starbucks!
* Starbucks is far away from the efficient frontier boundary

## Covariance between efficient portfolio returns

Later on, we will use the covariance between the two efficient portfolios.

```{r}
sigma.xy = as.numeric(t(x.vec)%*%sigma.mat%*%y.vec)
rho.xy = sigma.xy/(sig.px*sig.py)
sigma.xy
rho.xy
```

## Show efficient portfolios

```{r, echo=FALSE}
sd.vec = sqrt(diag(sigma.mat))
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)
points(sig.gmin, mu.gmin, pch=16, cex=2, col="green")
points(sig.px, mu.px, pch=16, cex=2, col="green")
points(sig.py, mu.py, pch=16, cex=2, col="green")      
text(sd.vec, mu.vec, labels=asset.names, pos=4, cex = cex.val)
text(sig.gmin, mu.gmin, labels="GLOBAL MIN", pos=2, cex = cex.val)        
text(sig.px, mu.px, labels="E1", pos=2, cex = cex.val) 
text(sig.py, mu.py, labels="E2", pos=2, cex = cex.val) 
text(0, r.f, labels=expression(r[f]), pos=2, cex = cex.val)

```

* Point E1 is the efficient portfolio with the same mean as Microsoft
* Point E2 is the efficient portfolio withe the same mean as Starbucks

## Find efficient portfolio from two efficient portfolios

Here we use the fact that any efficient portfolio is a convex combination of any two efficient portfolios: 

$$z = \alpha \times x + (1 - \alpha) \times y$$

Set $\alpha = 0.5$.

```{r}
a = 0.5
z.vec = a*x.vec + (1-a)*y.vec
z.vec
```

Compute the mean and volatility.

```{r}
sigma.xy = as.numeric(t(x.vec)%*%sigma.mat%*%y.vec)
mu.pz = as.numeric(crossprod(z.vec, mu.vec))
sig2.pz = as.numeric(t(z.vec)%*%sigma.mat%*%z.vec)
sig.pz = sqrt(sig2.pz)
mu.pz
sig.pz
```

Here the mean is half-way between the mean of Microsoft and the mean of Starbucks.

## Compute efficient portfolio with mean 0.05

Given a target mean value, $\mu_0 = 0.05$, you can solve for $\alpha$.

```{r}
a.05 = (0.05 - mu.py)/(mu.px - mu.py)
a.05
```

Given $\alpha=$ `r a.05` solve for $z$.

```{r}
z.05 = a.05*x.vec + (1 - a.05)*y.vec
z.05
```

Finally, compute the mean and volatility.

```{r}
mu.pz.05 = a.05*mu.px + (1-a.05)*mu.py
sig2.pz.05 = a.05^2 * sig2.px + (1-a.05)^2 * sig2.py + 2*a.05*(1-a.05)*sigma.xy          
sig.pz.05 = sqrt(sig2.pz.05)          
mu.pz.05
sig.pz.05
```


## Show efficient portfolios

```{r, echo=FALSE}
sd.vals = c(sig.px, sig.py, sig.pz, sig.pz.05)
mu.vals = c(mu.px, mu.py, mu.pz, mu.pz.05)
plot(sd.vals, mu.vals,  ylim=c(0, 0.06), xlim=c(0, 0.20), ylab=expression(mu[p]),
     xlab=expression(sigma[p]), pch=16, col="green", cex=2) 
points(sig.gmin, mu.gmin, pch=16, cex=2, col="green")     
points(sd.vec, mu.vec, pch=16, col="blue", cex=2)
points(sig.pz.05, mu.pz.05, pch=16, col="green", cex=2) 
text(sd.vals, mu.vals, labels=c("E1","E2","E3"), pos=2, cex = cex.val)
text(sig.gmin, mu.gmin, labels="GLOBAL MIN", pos=2, cex = cex.val)        
text(sd.vec, mu.vec, labels=asset.names, pos=4, cex = cex.val)
text(sig.pz.05, mu.pz.05, labels="E4", pos=2, cex = cex.val)  
text(0, r.f, labels=expression(r[f]), pos=2, cex = cex.val)
```

* Point E1 is the efficient portfolio with the same mean as Microsoft
* Point E2 is the efficient portfolio with the same mean as Starbucks
* Point E3 is the efficient portfolio with $\alpha = 0.5$
* Point E4 is the efficient portfolio with $\alpha=$ `r a.05`.

## Compute efficient frontier

Here we compute efficient portfolios as convex combinations of the the global minimum variance portfolio and the efficient portfolio with the same mean as Microsoft.

```{r}
a = seq(from=1, to=-1, by=-0.1)
n.a = length(a)
z.mat = matrix(0, n.a, 3)
colnames(z.mat) = names(mu.vec)
mu.z = rep(0, n.a)
sig2.z = rep(0, n.a)
sig.mx = t(m.vec)%*%sigma.mat%*%x.vec
for (i in 1:n.a) {
  z.mat[i, ] = a[i]*m.vec + (1-a[i])*x.vec
  mu.z[i] = a[i]*mu.gmin + (1-a[i])*mu.px
  sig2.z[i] = a[i]^2 * sig2.gmin + (1-a[i])^2 * sig2.px + 2*a[i]*(1-a[i])*sig.mx
}

```

## Show efficient frontier

```{r, echo=FALSE}
plot(sqrt(sig2.z), mu.z, type="b", ylim=c(0, 0.06), xlim=c(0, 0.16), 
     pch=16, col="green", cex = cex.val, ylab=expression(mu[p]), xlab=expression(sigma[p]))
points(sd.vec, mu.vec, pch=16, cex=2, lwd=2, col="blue")
points(sig.gmin, mu.gmin, pch=16, col="black", cex=2)
points(sig.px, mu.px, pch=16, col="black", cex=2)
text(sig.gmin, mu.gmin, labels="GLOBAL MIN", pos=2, cex = cex.val)
text(sd.vec, mu.vec, labels=asset.names, pos=4, cex = cex.val)
text(sig.px, mu.px, labels="E1", pos=2, cex = cex.val)
text(0, r.f, labels=expression(r[f]), pos=2, cex = cex.val)
```

## show weights in efficient portfolios

```{r}
chart.StackedBar(z.mat, xaxis.labels=round(sqrt(sig2.z),digits=3), 
                 xlab="Portfolio SD", ylab="Weights")
```

* The x-axis shows the volatility of each efficient portfolio
* The first portfolio is the global minimum variance portfolio
* Notice that Nordstrom eventually gets shorted in the efficient portfolios with high volatility and high expected return

## Show efficient frontier with random portfolios

```{r, echo=FALSE}
a = seq(from=1, to=-2, by=-0.1)
n.a = length(a)
z.mat = matrix(0, n.a, 3)
mu.z = rep(0, n.a)
sig2.z = rep(0, n.a)
sig.mx = t(m.vec)%*%sigma.mat%*%x.vec
for (i in 1:n.a) {
  z.mat[i, ] = a[i]*m.vec + (1-a[i])*x.vec
  mu.z[i] = a[i]*mu.gmin + (1-a[i])*mu.px
  sig2.z[i] = a[i]^2 * sig2.gmin + (1-a[i])^2 * sig2.px + 2*a[i]*(1-a[i])*sig.mx
}
plot(sqrt(sig2.z), mu.z, type="b", ylim=c(-0.03, 0.08), xlim=c(0, 0.4), 
     pch=16, col="green", cex = cex.val, ylab=expression(mu[p]), xlab=expression(sigma[p]))
points(sd.vec, mu.vec, pch=16, cex=2, col="blue")

for (i in 1:length(x.msft)) {
  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="black", cex=1.5)
}

text(sig.gmin, mu.gmin, labels="GLOBAL MIN", pos=2, col="black", cex = cex.val)
text(sd.vec, mu.vec, labels=asset.names, pos=4, col="blue", cex = cex.val)
text(0, r.f, labels=expression(r[f]), pos=2, cex = cex.val)

```

## Compute tangency (maximum Sharpe ratio) portfolio

Here we use the analytic matrix albegra formula for the tangency portfolio.

```{r}
rf = 0.005
sigma.inv.mat = solve(sigma.mat)
one.vec = rep(1, 3)
mu.minus.rf = mu.vec - rf*one.vec
top.mat = sigma.inv.mat%*%mu.minus.rf
bot.val = as.numeric(t(one.vec)%*%top.mat)
t.vec = top.mat[,1]/bot.val
t.vec
```

Compute mean and volatility of tangency portfolio

```{r}
mu.t = as.numeric(crossprod(t.vec, mu.vec))
sig2.t = as.numeric(t(t.vec)%*%sigma.mat%*%t.vec)
sig.t = sqrt(sig2.t)
mu.t
sig.t
```

## Show efficient portfolios when there is a risk free asset

Efficient portfolios are combinations of T-Bills and the tangency portfolio.

```{r, echo=FALSE}
x.t = seq(0, 2, by=0.1)
mu.pe = rf + x.t*(mu.t - rf)
sig.pe = x.t*sig.t 
slope.t = (mu.t - rf)/sig.t

# plot efficient portfolios
plot(sqrt(sig2.z), mu.z, type="b", ylim=c(0, 0.08), xlim=c(0, 0.17), 
     pch=16, col="blue", cex=1.15, ylab=expression(mu[p]), xlab=expression(sigma[p]))
abline(a=rf, b=slope.t, col="green", lwd=2)
points(sig.t, mu.t, pch=16, col="green", cex=1.15)
points(sd.vec, mu.vec, pch=16, cex=2, col="black")
text(sig.gmin, mu.gmin, labels="Global min", pos=4, cex=1.15)
text(sd.vec, mu.vec, labels=asset.names, pos=4, cex=1.15)
text(sig.t, mu.t, labels="tangency", pos=4,cex=1.15)
text(0, rf, labels=expression(r[f]), pos=2, cex=1.15)

```

## Show efficient portfolio weights

```{r, echo=FALSE}
t.mat = x.t %o% t.vec
e.mat = cbind(1-x.t, t.mat)
colnames(e.mat)[1] = "T-Bill"
chart.StackedBar(e.mat, xaxis.labels=round(sig.pe,digits=3), 
                 xlab="Portfolio SD", ylab="Weights")

```

## Find efficient portfolio with a target volatility

Every efficient portfolio is a combination of T-bills and the tangency portfolio. The volatility of such an efficient portfolio is: 

$$ \sigma_e = x_{tan} \times \sigma_{tan} $$

Given a target volatility, $\sigma_0 = 0.02$, you can solve for $x_{tan}$ and $x_f = 1 - x_{tan}$: 

```{r}
x.t.02 = 0.02/sig.t
x.t.02
1-x.t.02

```

The efficient portfolio weights in MSFT, NORD, SBUX are:

```{r}
x.t.02*t.vec
```

The mean and volatility of this efficient portfolio are:

```{r}
mu.t.02 = x.t.02*mu.t + (1-x.t.02)*rf
sig.t.02 = x.t.02*sig.t
mu.t.02
sig.t.02

```

## Find efficient portfolio with a target expected return

Every efficient portfolio is a combination of T-bills and the tangency portfolio. The mean of such an efficient portfolio is: 

$$ \mu_e = r_f + x_{tan} \times (\mu_{tan} - r_f) $$

Given a target mean, $\mu_0 = 0.07$, you can solve for $x_{tan}$ and $x_f = 1 - x_{tan}$: 

```{r}
x.t.07 = (0.07 - rf)/(mu.t - rf)
x.t.07
1-x.t.07
```

The efficient portfolio weights in MSFT, NORD, SBUX are:

```{r}
x.t.07*t.vec

```

The mean and volatility of this efficient portfolio are:

```{r}
mu.t.07 = x.t.07*mu.t + (1-x.t.07)*rf
sig.t.07 = x.t.07*sig.t
mu.t.07
sig.t.07

```

## Show efficient portfolios with target mean and target volatility

```{r, echo=FALSE}
plot(sqrt(sig2.z), mu.z, type="b", ylim=c(0, 0.08), xlim=c(0, 0.17), 
     pch=16, col="blue", cex=cex.val, ylab=expression(mu[p]), xlab=expression(sigma[p]))
abline(a=rf, b=slope.t, col="green", lwd=2)
points(sig.t, mu.t, pch=16, col="green", cex=2)
points(sig.t.02, mu.t.02, pch=16, col="green", cex=2)
points(sig.t.07, mu.t.07, pch=16, col="green", cex=2)
points(sd.vec, mu.vec, pch=16, col="black", cex=cex.val)
text(sig.gmin, mu.gmin, labels="Global min", pos=4, cex=cex.val)
text(sd.vec, mu.vec, labels=asset.names, pos=4, cex=cex.val)
text(sig.t, mu.t, labels="tangency", pos=4, cex=cex.val)
text(sig.t.02, mu.t.02, labels="E1", pos=3, cex=cex.val)
text(sig.t.07, mu.t.07, labels="E2", pos=3, cex=cex.val)
text(0, rf, labels=expression(r[f]), pos=2, cex=cex.val)

segments(0.02, 0, 0.02, mu.t.02, col="green", lty="dashed", lwd=2)
text(0.02, 0, labels="0.02", col="black", lwd=2, pos=3)
segments(0, mu.t.07, sig.t.07, mu.t.07, col="green", lty="dashed", lwd=2)
text(0, 0.07, labels="0.07", col="black", lwd=2, pos=4)

```

* Point E1 is the efficient portfolio with a target volatility $\sigma_0 = 0.02$.
* Point E2 is the efficient portfolio with a target mean $\mu_0 = 0.07$.
