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")
#install.packages("quadprog")
#install.packages("IntroCompFinR")
#install.packages("BLCOP")
#install.packages("Rdonlp2")
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)
library(quadprog)
library(fPortfolio)
## Loading required package: fBasics
## 
## Rmetrics Package fBasics
## Analysing Markets and calculating Basic Statistics
## Copyright (C) 2005-2014 Rmetrics Association Zurich
## Educational Software for Financial Engineering and Computational Science
## Rmetrics is free software and comes with ABSOLUTELY NO WARRANTY.
## https://www.rmetrics.org --- Mail to: info@rmetrics.org
## Loading required package: fAssets
## 
## Rmetrics Package fAssets
## Analysing and Modeling Financial Assets
## Copyright (C) 2005-2014 Rmetrics Association Zurich
## Educational Software for Financial Engineering and Computational Science
## Rmetrics is free software and comes with ABSOLUTELY NO WARRANTY.
## https://www.rmetrics.org --- Mail to: info@rmetrics.org
## 
## Rmetrics Package fPortfolio
## Portfolio Optimization
## Copyright (C) 2005-2014 Rmetrics Association Zurich
## Educational Software for Financial Engineering and Computational Science
## Rmetrics is free software and comes with ABSOLUTELY NO WARRANTY.
## https://www.rmetrics.org --- Mail to: info@rmetrics.org
## 
## Attaching package: 'fPortfolio'
## The following objects are masked from 'package:rrcov':
## 
##     getCov, getData
## The following object is masked from 'package:IntroCompFinR':
## 
##     getPortfolio
library(Rdonlp2)
## 
## Attaching package: 'Rdonlp2'
## The following object is masked from 'package:fPortfolio':
## 
##     donlp2NLP

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.

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
sd.vec = sqrt(diag(sigma.mat))
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

Random portfolios: no short sales

set.seed(123)
x.msft.0 = runif(400, min=-1.5, max=3)
x.nord.0 = runif(400, min=-1.5, max=3)
x.sbux.0 = 1 - x.msft.0 - x.nord.0

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)
## [1] 201

Risk-return characteristics: long-only random portfolios

Global minimum variance portfolio: no constraints on short sales

Use the IntroCompFinR package function globalMin.portfolio()

gmin.port = globalMin.portfolio(mu.vec, sigma.mat)
gmin.port
## Call:
## globalMin.portfolio(er = mu.vec, cov.mat = sigma.mat)
## 
## Portfolio expected return:     0.0249 
## Portfolio standard deviation:  0.0727 
## Portfolio weights:
##  MSFT  NORD  SBUX 
## 0.441 0.366 0.193

Use the quadprog function solve.QP() to compute global minimum variance portfolio

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

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))
D.mat
##        MSFT   NORD   SBUX
## MSFT 0.0200 0.0036 0.0022
## NORD 0.0036 0.0218 0.0052
## SBUX 0.0022 0.0052 0.0398
d.vec
## [1] 0 0 0
t(A.mat)
##      [,1] [,2] [,3]
## [1,]   -1   -1   -1
## [2,]    1    0    0
## [3,]    0    1    0
## [4,]    0    0    1
b.vec
## [1] -1  0  0  0

Use the quadprog function solve.QP() to compute global minimum variance portfolio

The function solve.QP() takes the restrictions matrices as inputs

args(solve.QP)
## function (Dmat, dvec, Amat, bvec, meq = 0, factorized = FALSE) 
## NULL

solve.QP() returns a list containing information about the optimization

qp.out = solve.QP(Dmat=D.mat, dvec=d.vec,
                  Amat=A.mat, bvec=b.vec, meq=1)
class(qp.out)
## [1] "list"

The components of the returned list are

names(qp.out)
## [1] "solution"               "value"                 
## [3] "unconstrained.solution" "iterations"            
## [5] "Lagrangian"             "iact"

The portfolio weights are in the solution component

qp.out$solution
## [1] 0.441 0.366 0.193

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

sum(qp.out$solution)
## [1] 1

The minimized value of the objective function (portfolio variance) is in the value component

qp.out$value
## [1] 0.00528

Use the quadprog function solve.QP() to compute global minimum variance portfolio

Compute the mean and volatility of the minimum variance portfolio

# Expected return
er.gmin.ns = as.numeric(crossprod(qp.out$solution, mu.vec))
# Volatility
sd.gmin.ns = sqrt(qp.out$value)
er.gmin.ns
## [1] 0.0249
sd.gmin.ns
## [1] 0.0727

Use the IntroCompFinR function globalMin.Portfolio() to compute the minimum variance portfolio

Use the optional argument shorts=TRUE to impose no-short sales restrictions

args(globalMin.portfolio)
## function (er, cov.mat, shorts = TRUE) 
## NULL

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

gmin.port = globalMin.portfolio(mu.vec, sigma.mat,
                                shorts=FALSE)
gmin.port
## Call:
## globalMin.portfolio(er = mu.vec, cov.mat = sigma.mat, shorts = FALSE)
## 
## Portfolio expected return:     0.0249 
## Portfolio standard deviation:  0.0727 
## Portfolio weights:
##  MSFT  NORD  SBUX 
## 0.441 0.366 0.193

Compute minimum variance portfolio with same mean as Microsoft

First find the minimum variance portfolio allowing short sales

eMsft.port = efficient.portfolio(mu.vec, sigma.mat, target.return = mu.vec["MSFT"])
eMsft.port
## Call:
## efficient.portfolio(er = mu.vec, cov.mat = sigma.mat, target.return = mu.vec["MSFT"])
## 
## Portfolio expected return:     0.0427 
## Portfolio standard deviation:  0.0917 
## Portfolio weights:
##    MSFT    NORD    SBUX 
##  0.8275 -0.0907  0.2633

Compute minimum variance portfolio with same mean as Microsoft

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

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))
D.mat
##        MSFT   NORD   SBUX
## MSFT 0.0200 0.0036 0.0022
## NORD 0.0036 0.0218 0.0052
## SBUX 0.0022 0.0052 0.0398
d.vec
## [1] 0 0 0
A.mat
##      mu.vec        
## MSFT 0.0427 1 1 0 0
## NORD 0.0015 1 0 1 0
## SBUX 0.0285 1 0 0 1
b.vec
##   MSFT                             
## 0.0427 1.0000 0.0000 0.0000 0.0000

Compute minimum variance portfolio with same mean as Microsoft

Next we use solve.QP() to find the solution

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)
## MSFT NORD SBUX 
##    1    0    0

Compute minimum variance portfolio with same mean as Microsoft

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

efficient.portfolio(mu.vec, sigma.mat, target.return = mu.vec["MSFT"], shorts = FALSE)
## Call:
## efficient.portfolio(er = mu.vec, cov.mat = sigma.mat, target.return = mu.vec["MSFT"], 
##     shorts = FALSE)
## 
## Portfolio expected return:     0.0427 
## Portfolio standard deviation:  0.1 
## Portfolio weights:
## MSFT NORD SBUX 
##    1    0    0

Infeasible no-short sale efficient portfolio

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

# 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!

Compute efficient frontier allowing for short sales

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

ef = efficient.frontier(mu.vec, sigma.mat, alpha.min=0, 
                         alpha.max=1, nport=10)
ef$weights
##          MSFT    NORD  SBUX
## port 1  0.827 -0.0907 0.263
## port 2  0.785 -0.0400 0.256
## port 3  0.742  0.0107 0.248
## port 4  0.699  0.0614 0.240
## port 5  0.656  0.1121 0.232
## port 6  0.613  0.1628 0.224
## port 7  0.570  0.2135 0.217
## port 8  0.527  0.2642 0.209
## port 9  0.484  0.3149 0.201
## port 10 0.441  0.3656 0.193

Compute efficient frontier allowing for short sales

Compute efficient frontier without short sales

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

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

Compute efficient frontier without short sales

ef.ns = efficient.frontier(mu.vec, sigma.mat, alpha.min=0, 
                         alpha.max=1, nport=10, shorts=FALSE)
ef.ns$weights
##          MSFT   NORD  SBUX
## port 1  0.441 0.3656 0.193
## port 2  0.484 0.3149 0.201
## port 3  0.527 0.2642 0.209
## port 4  0.570 0.2135 0.217
## port 5  0.613 0.1628 0.224
## port 6  0.656 0.1121 0.232
## port 7  0.699 0.0614 0.240
## port 8  0.742 0.0107 0.248
## port 9  0.861 0.0000 0.139
## port 10 1.000 0.0000 0.000

Compare efficient frontiers

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 1", col="blue", pos=3)
text(ef$sd[2], ef$er[2], labels="port 2", col="blue", pos=3)
text(ef.ns$sd[10], ef.ns$er[10], labels="port 1", col="red", pos=1)
text(ef.ns$sd[9], ef.ns$er[9], labels="port 2", col="red", pos=1)

library(xlsx)
## Loading required package: rJava
## Loading required package: xlsxjars
temp <- tempfile()
download.file("http://freakonometrics.free.fr/oil.xls",temp)
oil=read.xlsx(temp,sheetName="DATA",dec=",")
head(oil)
##         Date    WTI  brent  Dubai   Maya NA. dec
## 1 1997-01-10  2.737  2.255  3.367   1.54  NA   ,
## 2 1997-01-17 -3.403 -6.014 -3.825  -4.11  NA   ,
## 3 1997-01-24 -4.095 -1.431 -6.638  -4.62  NA   ,
## 4 1997-01-31 -0.658  0.349  0.733  -1.51  NA   ,
## 5 1997-02-07 -3.143 -1.978 -0.733  -1.88  NA   ,
## 6 1997-02-14 -5.603 -7.845 -7.637 -11.05  NA   ,
Time=as.Date(oil$Date,"%Y-%m-%d")
plot(Time,oil[,3],type="l",ylab="Brent, weekly log returns",ylim=range(oil[,3:5]))