Set options and load packages

This course is inspired by notes by Eric Zivot’s lecture notes and Bernhard Pfaff’s Financial Risk Modelling and Portfolio Optimization with R. Slides are available from my blog. See https://freakonometrics.hypotheses.org/ for more details.

options(digits=3, width=70)
install.packages("IntroCompFinR", repos="http://R-Forge.R-project.org")
install.packages("PerformanceAnalytics")
install.packages("zoo")
install.packages("zoo")
install.packages("rrcov")
install.packages("FRAPO")
Sys.setenv(TZ="UTC")
options(digits=3, width=70)
cex.val = 2
library(IntroCompFinR)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(PerformanceAnalytics)
## 
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
library(zoo)
library(zoo)
library(rrcov)
## Loading required package: robustbase
## Scalable Robust Estimators with High Breakdown Point (version 1.4-3)
library(FRAPO)
## Loading required package: cccp
## Loading required package: Rglpk
## Loading required package: slam
## Using the GLPK callable library version 4.61
## Loading required package: timeSeries
## Loading required package: timeDate
## 
## Attaching package: 'timeDate'
## The following objects are masked from 'package:PerformanceAnalytics':
## 
##     kurtosis, skewness
## 
## Attaching package: 'timeSeries'
## The following object is masked from 'package:zoo':
## 
##     time<-
## Financial Risk Modelling and Portfolio Optimisation with R (version 0.4-1)

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.

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
##   MSFT   NORD   SBUX 
## 0.0427 0.0015 0.0285
cov2cor(sigma.mat)
##       MSFT  NORD  SBUX
## MSFT 1.000 0.172 0.078
## NORD 0.172 1.000 0.177
## SBUX 0.078 0.177 1.000

Risk return characteristics

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=1, cex.lab=1.75)     
text(sd.vec, mu.vec, labels=asset.names, pos=4, cex = 1.25)
text(0, r.f, labels=expression(r[f]), pos=2, cex = 1.25)

Example portfolio: equally weighted

# 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)
# show mean and volatility
mu.p.x
##        [,1]
## [1,] 0.0242
sig.p.x
##        [,1]
## [1,] 0.0759

Example portfolio: long-short

# 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)
# show mean and volatility
mu.p.y
##        [,1]
## [1,] 0.0291
sig.p.y
##        [,1]
## [1,] 0.0966

Covariance and correlation between example portfolio returns

# covariance
sig.xy = t(x.vec)%*%sigma.mat%*%y.vec
sig.xy
##         [,1]
## [1,] 0.00391
# correlation
rho.xy = sig.xy/(sig.p.x*sig.p.y)
rho.xy
##       [,1]
## [1,] 0.533

Risk-return characteristics: example portfolios

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 and plot them on the volatility-mean plot

set.seed(1)
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))
##      x.msft x.nord x.sbux
## [1,] -0.703  0.464  1.239
## [2,] -0.384 -0.440  1.824
## [3,]  0.219 -0.689  1.471
## [4,]  1.225  1.478 -1.703
## [5,] -0.895  0.400  1.494
## [6,]  1.195 -0.860  0.665

Risk-return characteristics: random portfolios

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)
}

Compute global minimum variance portfolio: Method 1

Use \(z_m = A^{-1}_m b\).

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
##  MSFT  NORD  SBUX 
## 0.441 0.366 0.193

Mean and volatility of minimum variance portfolio

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
## [1] 0.0249
sig.gmin
## [1] 0.0727

Compute global minimum variance portfolio: Method 2

Use analytic formula for minimum variance portfolio

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]
##  MSFT  NORD  SBUX 
## 0.441 0.366 0.193

Show minimum variance portfolio

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)
}

Compute efficient portfolio with the same mean as Microsoft

Use matrix algebra formula to compute efficient portfolio.

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
##    MSFT    NORD    SBUX 
##  0.8275 -0.0907  0.2633

Compute mean and volatility of efficient portfolio.

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
## [1] 0.0427
sig.px
## [1] 0.0917

Compare with mean and volatility of Microsoft.

mu.vec["MSFT"]
##   MSFT 
## 0.0427
sd.vec["MSFT"]
## MSFT 
##  0.1

Compute efficient portfolio with the same mean as Starbucks

# 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
##  MSFT  NORD  SBUX 
## 0.519 0.273 0.207
# 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
## [1] 0.0285
sig.py
## [1] 0.0736
# compare with Starbucks
mu.vec["SBUX"]
##   SBUX 
## 0.0285
sd.vec["SBUX"]
##  SBUX 
## 0.141

Covariance between efficient portfolio returns

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

sigma.xy = as.numeric(t(x.vec)%*%sigma.mat%*%y.vec)
rho.xy = sigma.xy/(sig.px*sig.py)
sigma.xy
## [1] 0.00591
rho.xy
## [1] 0.877

Show efficient portfolios

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)

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

a = 0.5
z.vec = a*x.vec + (1-a)*y.vec
z.vec
##   MSFT   NORD   SBUX 
## 0.6734 0.0912 0.2354

Compute the mean and volatility.

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
## [1] 0.0356
sig.pz
## [1] 0.0801

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

a.05 = (0.05 - mu.py)/(mu.px - mu.py)
a.05
## [1] 1.51

Given \(\alpha=\) 1.514 solve for \(z\).

z.05 = a.05*x.vec + (1 - a.05)*y.vec
z.05
##   MSFT   NORD   SBUX 
##  0.986 -0.278  0.292

Finally, compute the mean and volatility.

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
## [1] 0.05
sig.pz.05
## [1] 0.107

Show efficient portfolios

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)