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)

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.

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

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

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

Show efficient frontier with random portfolios

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.

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
##   MSFT   NORD   SBUX 
##  1.027 -0.326  0.299

Compute mean and volatility of tangency portfolio

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
## [1] 0.0519
sig.t
## [1] 0.112

Show efficient portfolios when there is a risk free asset

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

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(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

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

x.t.02 = 0.02/sig.t
x.t.02
## [1] 0.179
1-x.t.02
## [1] 0.821

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

x.t.02*t.vec
##    MSFT    NORD    SBUX 
##  0.1840 -0.0585  0.0537

The mean and volatility of this efficient portfolio are:

mu.t.02 = x.t.02*mu.t + (1-x.t.02)*rf
sig.t.02 = x.t.02*sig.t
mu.t.02
## [1] 0.0134
sig.t.02
## [1] 0.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}\):

x.t.07 = (0.07 - rf)/(mu.t - rf)
x.t.07
## [1] 1.39
1-x.t.07
## [1] -0.386

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

x.t.07*t.vec
##   MSFT   NORD   SBUX 
##  1.423 -0.452  0.415

The mean and volatility of this efficient portfolio are:

mu.t.07 = x.t.07*mu.t + (1-x.t.07)*rf
sig.t.07 = x.t.07*sig.t
mu.t.07
## [1] 0.07
sig.t.07
## [1] 0.155

Show efficient portfolios with target mean and target volatility

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)

Real Time Series

data( StockIndex )
pzoo = zoo ( StockIndex , order.by = rownames ( StockIndex ) )
rzoo = ( pzoo / lag ( pzoo , k = -1) - 1 ) * 100

xx

Moments <- function ( x , method = c ( "CovClassic" , "CovMcd" , "CovMest" , "CovMMest" , "CovMve" , "CovOgk" , "CovSde" , "CovSest" ) , ... ) {
method <- match.arg ( method )
ans <- do.call ( method , list ( x = x , ... ) )
return ( getCov ( ans ) )}
Moments(as.matrix(rzoo),"CovClassic")
##         SP500 N225 FTSE100 CAC40 GDAX  HSI
## SP500    17.8 12.7    13.8  17.8 19.5 18.9
## N225     12.7 36.6    10.8  15.0 16.2 16.7
## FTSE100  13.8 10.8    17.3  18.8 19.4 19.1
## CAC40    17.8 15.0    18.8  30.9 29.9 22.8
## GDAX     19.5 16.2    19.4  29.9 38.0 26.1
## HSI      18.9 16.7    19.1  22.8 26.1 58.1
Moments(as.matrix(rzoo),"CovMcd")
##         SP500 N225 FTSE100 CAC40 GDAX  HSI
## SP500    18.0 14.8    13.8  17.8 20.2 17.4
## N225     14.8 41.1    12.9  18.4 18.9 15.7
## FTSE100  13.8 12.9    16.1  17.6 19.3 14.9
## CAC40    17.8 18.4    17.6  28.3 27.6 17.0
## GDAX     20.2 18.9    19.3  27.6 34.5 21.7
## HSI      17.4 15.7    14.9  17.0 21.7 44.1
Moments(as.matrix(rzoo),"CovMest")
##         SP500 N225 FTSE100 CAC40 GDAX  HSI
## SP500    16.9 13.1    12.8  16.6 18.5 16.5
## N225     13.1 39.7    11.2  15.9 16.9 16.5
## FTSE100  12.8 11.2    14.9  16.2 17.3 13.9
## CAC40    16.6 15.9    16.2  26.7 25.4 16.1
## GDAX     18.5 16.9    17.3  25.4 31.6 20.9
## HSI      16.5 16.5    13.9  16.1 20.9 43.1
Moments(as.matrix(rzoo),"CovMMest")
##         SP500 N225 FTSE100 CAC40 GDAX  HSI
## SP500    16.1 12.0    12.3  16.0 17.7 16.3
## N225     12.0 35.6    10.3  14.6 15.4 14.5
## FTSE100  12.3 10.3    15.1  16.4 17.4 15.4
## CAC40    16.0 14.6    16.4  27.6 26.4 18.1
## GDAX     17.7 15.4    17.4  26.4 33.0 21.8
## HSI      16.3 14.5    15.4  18.1 21.8 46.6
Moments(as.matrix(rzoo),"CovOgk")
##         SP500  N225 FTSE100 CAC40 GDAX  HSI
## SP500   12.76  9.87    9.38  12.3 13.6 12.3
## N225     9.87 27.53    8.25  11.9 12.6 11.4
## FTSE100  9.38  8.25   10.92  11.8 12.4 10.4
## CAC40   12.32 11.91   11.81  19.7 18.6 12.3
## GDAX    13.61 12.63   12.42  18.6 23.2 15.0
## HSI     12.29 11.45   10.37  12.3 15.0 33.0
Moments(as.matrix(rzoo),"CovSde")
##         SP500 N225 FTSE100 CAC40 GDAX  HSI
## SP500    16.9 14.4    12.1  16.1 18.5 18.4
## N225     14.4 39.7    12.3  18.3 19.4 16.7
## FTSE100  12.1 12.3    13.9  14.9 16.6 15.3
## CAC40    16.1 18.3    14.9  25.5 25.2 18.2
## GDAX     18.5 19.4    16.6  25.2 31.8 22.7
## HSI      18.4 16.7    15.3  18.2 22.7 47.0

Computations Based on Those Matrices

download.file(url="http://freakonometrics.free.fr/portfolio.r",destfile = "portfolio.r")
source("portfolio.r")

Compute various estimates of variance matrices, for those six indices

er=apply(rzoo,2,mean)
covmat=Moments(as.matrix(rzoo),"CovClassic")
globalMin.portfolio(er, covmat)
## Call:
## globalMin.portfolio(er = er, cov.mat = covmat)
## 
## Portfolio expected return:     0.306 
## Portfolio standard deviation:  3.71 
## Portfolio weights:
##   SP500    N225 FTSE100   CAC40    GDAX     HSI 
##  0.5220  0.1531  0.7196 -0.1461 -0.1675 -0.0811
covmat=Moments(as.matrix(rzoo),"CovMcd")
globalMin.portfolio(er, covmat)
## Call:
## globalMin.portfolio(er = er, cov.mat = covmat)
## 
## Portfolio expected return:     0.309 
## Portfolio standard deviation:  3.65 
## Portfolio weights:
##   SP500    N225 FTSE100   CAC40    GDAX     HSI 
##  0.5309  0.0754  0.8365 -0.0197 -0.4176 -0.0055
covmat=Moments(as.matrix(rzoo),"CovMest")
globalMin.portfolio(er, covmat)
## Call:
## globalMin.portfolio(er = er, cov.mat = covmat)
## 
## Portfolio expected return:     0.32 
## Portfolio standard deviation:  3.59 
## Portfolio weights:
##   SP500    N225 FTSE100   CAC40    GDAX     HSI 
##  0.4827  0.0976  0.8168 -0.0671 -0.3193 -0.0107
covmat=Moments(as.matrix(rzoo),"CovMMest")
globalMin.portfolio(er, covmat)
## Call:
## globalMin.portfolio(er = er, cov.mat = covmat)
## 
## Portfolio expected return:     0.316 
## Portfolio standard deviation:  3.52 
## Portfolio weights:
##   SP500    N225 FTSE100   CAC40    GDAX     HSI 
##  0.5225  0.1216  0.7513 -0.1002 -0.2448 -0.0504
covmat=Moments(as.matrix(rzoo),"CovOgk")
globalMin.portfolio(er, covmat)
## Call:
## globalMin.portfolio(er = er, cov.mat = covmat)
## 
## Portfolio expected return:     0.322 
## Portfolio standard deviation:  3.08 
## Portfolio weights:
##   SP500    N225 FTSE100   CAC40    GDAX     HSI 
##  0.4441  0.1077  0.8278 -0.0963 -0.2637 -0.0196
covmat=Moments(as.matrix(rzoo),"CovSde")
globalMin.portfolio(er, covmat)
## Call:
## globalMin.portfolio(er = er, cov.mat = covmat)
## 
## Portfolio expected return:     0.29 
## Portfolio standard deviation:  3.46 
## Portfolio weights:
##   SP500    N225 FTSE100   CAC40    GDAX     HSI 
##  0.5013  0.0505  0.8830  0.0393 -0.3955 -0.0785
target.return <- 1
e.port.msft <- efficient.portfolio(er, covmat, target.return)
e.port.msft
## Call:
## efficient.portfolio(er = er, cov.mat = covmat, target.return = target.return)
## 
## Portfolio expected return:     1 
## Portfolio standard deviation:  4.84 
## Portfolio weights:
##   SP500    N225 FTSE100   CAC40    GDAX     HSI 
##   0.725  -0.380   0.476  -0.322   0.344   0.157

We can also derive the tangency portfolio, given some risk free asset at some given rate

target.return <- 1
tan.port <- tangency.portfolio(er, covmat, .06)
tan.port
## Call:
## tangency.portfolio(er = er, cov.mat = covmat, risk.free = 0.06)
## 
## Portfolio expected return:     2.57 
## Portfolio standard deviation:  11.4 
## Portfolio weights:
##   SP500    N225 FTSE100   CAC40    GDAX     HSI 
##   1.219  -1.330  -0.423  -1.119   1.977   0.677

or the efficient frontier

ef <- efficient.frontier(er, covmat, alpha.min=-2, alpha.max=1.5, nport=20)
plot(ef)
points(apply(rzoo,2,sd),apply(rzoo,2,mean),pch=19,col="blue")