Optimal Portfolio #2

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

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

• The argument meq = 2 specifies 2 equality constraints. It tells solve.QP() 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 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!

• Notice that solve.QP() throws an error saying there is no compatible 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

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

Compute efficient frontier allowing for short sales

• 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

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

• Here, the no-short sales efficient frontier lies “inside” and “to the right” of the short sales frontier

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

dat=oil[,c(3,4,5)]

ARMA

fit1=arima(x=dat[,1],order=c(2,0,1))
fit2=arima(x=dat[,2],order=c(1,0,1))
fit3=arima(x=dat[,3],order=c(1,0,1))
dat_arma <- cbind(residuals(fit1),residuals(fit2),residuals(fit3))
m <- apply(dat_arma, 2, mean)
v <- apply(dat_arma, 2, var)
dat_arma_std <- t((t(dat_arma)-m)/sqrt(v))

library(fGarch)
fit1=garchFit(formula = ~ arma(2,1)+garch(1, 1),data=dat[,1],cond.dist ="std",trace=FALSE)
fit2=garchFit(formula = ~ arma(1,1)+garch(1, 1),data=dat[,2],cond.dist ="std",trace=FALSE)
fit3=garchFit(formula = ~ arma(1,1)+garch(1, 1),data=dat[,3],cond.dist ="std",trace=FALSE)
dat_res<- cbind(residuals(fit1),residuals(fit2),residuals(fit3))
m_res <- apply(dat_res, 2, mean)
v_res <- apply(dat_res, 2, var)
dat_res_std =cbind((dat_res[,1]-m_res[1])/sqrt(v_res[1]),(dat_res[,2]-m_res[2])/sqrt(v_res[2]),(dat_res[,3]-m_res[3])/sqrt(v_res[3]))

library(MTS)
ewma=EWMAvol(dat_res_std, lambda = 0.96)
emwa_series_vol=function(i=1){
 plot(Time,ewma$Sigma.t[,i],type="l",
 col="white",ylim=c(-5,3))
 lines(Time,dat_res_std[,i],col="grey")
 j=diag(matrix(1:9,3,3))[i]
 lines(Time,ewma$Sigma.t[,j],lwd=1.5)
}
emwa_series_vol(1)

emwa_series_vol(2)

emwa_series_vol(3)

emwa_series_cor=function(i=1,j=2){
 if((min(i,j)==1)&(max(i,j)==2)){
   a=1;     b=5;     ab=2}
 if((min(i,j)==1)&(max(i,j)==3)){
   a=1;     b=9;     ab=3}
 if((min(i,j)==2)&(max(i,j)==3)){
   a=5;     b=9;     ab=6}
r=ewma$Sigma.t[,ab]/sqrt(ewma$Sigma.t[,a]*
ewma$Sigma.t[,b])
plot(Time,r,type="l",ylim=c(0,1))
}
emwa_series_cor(1,2)

emwa_series_cor(1,3)

emwa_series_cor(2,3)

library(MTS)
bekk=BEKK11(dat_arma)

## Initial estimates:  -0.000118 -0.000681 -0.00139 4.56 3.07 4.26 2.67 0.999 2.73 0.1 0.02 0.02 0.02 0.1 0.02 0.02 0.02 0.1 0.8 0.02 0.02 0.02 0.8 0.02 0.02 0.02 0.8 
## Lower limits:  -0.00118 -0.00681 -0.0139 0.912 0.615 0.851 0.533 0.2 0.545 1e-06 -0.5 -0.5 -0.5 1e-06 -0.5 -0.5 -0.5 1e-06 1e-06 -0.5 -0.5 -0.5 1e-06 -0.5 -0.5 -0.5 1e-06 
## Upper limits:  0.00118 0.00681 0.0139 5.02 3.38 4.68 2.93 1.1 3 1 0.5 0.5 0.5 1 0.5 0.5 0.5 1 1 0.5 0.5 0.5 1 0.5 0.5 0.5 1

## Warning in sqrt(diag(solve(Hessian))): production de NaN

## 
## Coefficient(s):
##                         Estimate   Std. Error  t value   Pr(>|t|)    
## mu1.residuals(fit1)  0.001184943  0.146640656  0.00808 0.99355269    
## mu2.residuals(fit2)  0.006811456  0.139468358  0.04884 0.96104783    
## mu3.residuals(fit3) -0.013904020  0.161408891 -0.08614 0.93135386    
## A011                 4.420544400  0.300451502 14.71300 < 2.22e-16 ***
## A021                 3.226348218  0.700553711  4.60543 4.1162e-06 ***
## A031                 4.439810974  0.434352154 10.22169 < 2.22e-16 ***
## A022                 0.533417942  2.195688225  0.24294 0.80805278    
## A032                 0.296108081  1.230932016  0.24056 0.80989926    
## A033                 0.545385862  0.134268540  4.06190 4.8674e-05 ***
## A11                  0.037699169  0.121109724  0.31128 0.75558693    
## A21                 -0.488770235  0.100894066 -4.84439 1.2700e-06 ***
## A31                 -0.440918277  0.125514239 -3.51289 0.00044325 ***
## A12                  0.253097413  0.109618095  2.30890 0.02094903 *  
## A22                  0.523246131  0.094432266  5.54097 3.0080e-08 ***
## A32                  0.272907075  0.116204861  2.34850 0.01884922 *  
## A13                 -0.054333597  0.120784788 -0.44984 0.65282720    
## A23                 -0.059993023  0.099765199 -0.60134 0.54761209    
## A33                  0.345443272  0.128314036  2.69217 0.00709886 ** 
## B11                  0.000001000  0.000238529  0.00419 0.99665499    
## B21                  0.500000000  0.204920750  2.43997 0.01468858 *  
## B31                 -0.291410720  0.172109283 -1.69317 0.09042262 .  
## B12                  0.024871918  0.362416448  0.06863 0.94528573    
## B22                  0.000001000           NA       NA         NA    
## B32                 -0.425339912  0.219554189 -1.93729 0.05271003 .  
## B13                 -0.127741634  0.241423728 -0.52912 0.59672361    
## B23                 -0.146411810  0.112278047 -1.30401 0.19222983    
## B33                  0.736846625  0.198298739  3.71584 0.00020253 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

bekk_series_vol=function(i=1){
   plot(Time,bekk$Sigma.t[,1],type="l",
        ylab=names(dat)[i],col="white",ylim=c(0,80))
   lines(Time,dat_arma[,i]+40,col="grey")
 j=diag(matrix(1:9,3,3))[i]
   lines(Time,bekk$Sigma.t[,j],lwd=1.5)
}
bekk_series_vol(1)

bekk_series_vol(2)

bekk_series_vol(3)

bekk_series_cor=function(i=1,j=2){
 if((min(i,j)==1)&(max(i,j)==2)){
   a=1;     b=5;     ab=2}
 if((min(i,j)==1)&(max(i,j)==3)){
   a=1;     b=9;     ab=3}
 if((min(i,j)==2)&(max(i,j)==3)){
   a=5;     b=9;     ab=6}
 r=bekk$Sigma.t[,ab]/sqrt(bekk$Sigma.t[,a]*bekk$Sigma.t[,b])
 plot(Time,r,type="l",
      ylim=c(0,1))
}
bekk_series_cor(1,2)

bekk_series_cor(1,3)

bekk_series_cor(2,3)

time_varying_correl_2=function(i=1,j=2,nom_arg="Pearson"){
   uv=dat_arma[,c(i,j)]
   n=nrow(uv)
   correlation=function(t){
       sub_uv=uv[t+(-26):26,]
       cor(sub_uv,method = tolower(nom_arg))[1,2]
   }
   ic_correlation=function(t){
       b=1000
       sub_uv=uv[t+(-26):26,]
       cr=rep(NA,b)
       for(sim in 1:b){
           sub_uv_b=sub_uv[sample(1:53,size=53,replace=TRUE),]+matrix(rnorm(106),53,2)/10000
           cr[sim]=cor(sub_uv_b,method = tolower(nom_arg))[1,2]}
       quantile(cr,c(.05,.95))
   }
   nom=paste("Corr?lation de ",nom_arg, " ",names(dat)[i]," vs. ",names(dat)[j],sep="")
   C=sapply(26:(n-26),correlation)
   IC_C=sapply(27:(n-27),ic_correlation)
   plot(Time[26:(n-26)],C,type="l",xlab="Time",ylab=nom,ylim=c(.3,1),col="white")
   polygon(c(Time[27:(n-27)],rev(Time[27:(n-27)])),c(IC_C[1,],rev(IC_C[2,])),col="light blue",border=NA)  
   lines(Time[26:(n-26)],C)  
   abline(h=cor(uv,method = tolower(nom_arg))[1,2],lty=2,col="red")
}

time_varying_correl_2(1,2)

time_varying_correl_2(1,3)

time_varying_correl_2(2,3)

library(MTS)
m2=dccFit(dat_res_std)

## Estimates:  0.92 0.0276 5.14 
## st.errors:  0.03 0.00947 0.428 
## t-values:   30.6 2.91 12

m3=dccFit(dat_res_std,type="Engle")

## Warning in sqrt(diag(Hi)): production de NaN

## Estimates:  0.92 0.0254 5.49 
## st.errors:  NaN 0.00778 0.467 
## t-values:   NaN 3.27 11.8

R2=m2$rho.t
plot(Time,R2[,2],type="l",ylim=c(0,1))

plot(Time,R2[,3],type="l",ylim=c(0,1))

plot(Time,R2[,4],type="l",ylim=c(0,1))

R3=m3$rho.t
plot(Time,R3[,2],type="l",ylim=c(0,1))

plot(Time,R3[,3],type="l",ylim=c(0,1))

plot(Time,R3[,4],type="l",ylim=c(0,1))

Back on the Stock Dataset

library(zoo)
data( StockIndex )
pzoo = zoo ( StockIndex , order.by = rownames ( StockIndex ) )
rzoo = ( pzoo / lag ( pzoo , k = -1) - 1 ) 
dat=rzoo
fit1=arima(x=dat[,1],order=c(1,0,1))

## Warning in as.ts.zoo(x): 'x' does not have an underlying regularity

fit2=arima(x=dat[,2],order=c(1,0,1))

## Warning in as.ts.zoo(x): 'x' does not have an underlying regularity

fit3=arima(x=dat[,3],order=c(1,0,1))

## Warning in as.ts.zoo(x): 'x' does not have an underlying regularity

fit4=arima(x=dat[,4],order=c(1,0,1))

## Warning in as.ts.zoo(x): 'x' does not have an underlying regularity

fit5=arima(x=dat[,5],order=c(1,0,1))

## Warning in as.ts.zoo(x): 'x' does not have an underlying regularity

fit6=arima(x=dat[,6],order=c(1,0,1))

## Warning in as.ts.zoo(x): 'x' does not have an underlying regularity

library(fGarch)
fit1=garchFit(formula = ~ arma(1,1)+garch(1, 1),data=dat[,1],cond.dist ="std",trace=FALSE)

## Warning in sqrt(diag(fit$cvar)): production de NaN

fit2=garchFit(formula = ~ arma(1,1)+garch(1, 1),data=dat[,2],cond.dist ="std",trace=FALSE)
fit3=garchFit(formula = ~ arma(1,1)+garch(1, 1),data=dat[,3],cond.dist ="std",trace=FALSE)
fit4=garchFit(formula = ~ arma(1,1)+garch(1, 1),data=dat[,1],cond.dist ="std",trace=FALSE)

## Warning in sqrt(diag(fit$cvar)): production de NaN

fit5=garchFit(formula = ~ arma(1,1)+garch(1, 1),data=dat[,2],cond.dist ="std",trace=FALSE)
fit6=garchFit(formula = ~ arma(1,1)+garch(1, 1),data=dat[,3],cond.dist ="std",trace=FALSE)
dat_res<- cbind(residuals(fit1),residuals(fit2),residuals(fit3),residuals(fit4),residuals(fit5),residuals(fit6))
m_res <- apply(dat_res, 2, mean)
v_res <- apply(dat_res, 2, var)
dat_res_std =cbind((dat_res[,1]-m_res[1])/sqrt(v_res[1]),(dat_res[,2]-m_res[2])/sqrt(v_res[2]),(dat_res[,3]-m_res[3])/sqrt(v_res[3]),(dat_res[,4]-m_res[4])/sqrt(v_res[4]),(dat_res[,5]-m_res[5])/sqrt(v_res[5]),(dat_res[,6]-m_res[6])/sqrt(v_res[6]))

Time=time(dat)
library(MTS)
ewma=EWMAvol(dat_res_std, lambda = 0.96)
emwa_series_vol=function(i=1){
 plot(ts(ewma$Sigma.t[,i],start=as.Date("1991/08/30","%Y/%m/%d"),frequency=12),type="l",
 col="white",ylim=c(-4,4),ylab="")
 lines(ts(dat_res_std[,i],start=as.Date("1991/08/30","%Y/%m/%d"),frequency=12),col="grey")
  j=diag(matrix(1:36,6,6))[i]
 lines(ts(ewma$Sigma.t[,j],start=as.Date("1991/08/30","%Y/%m/%d"),frequency=12)-4,lwd=1.5)
}
emwa_series_vol(1)

emwa_series_vol(2)

emwa_series_vol(3)

emwa_series_vol(4)

emwa_series_vol(5)

emwa_series_vol(6)

Time varying means and variances (and its impact on minimal variance portfolio)

k=nrow(rzoo)-60
P=matrix(NA,k,6)
for(j in 1:k){
covmat=var(as.matrix(rzoo)[j-1+1:60,])
er=apply(as.matrix(rzoo[j-1+1:60,]),2,mean)
P[j,]=globalMin.portfolio (er , covmat , shorts = FALSE )$weights
}
 library(RColorBrewer)
 darkcols <- brewer.pal(8, "Dark2")
 add_legend <- function(...) {
   opar <- par(fig=c(0, 1, 0, 1), oma=c(0, 0, 0, 0), 
               mar=c(0, 0, 0, 0), new=TRUE)
   on.exit(par(opar))
   plot(0, 0, type='n', bty='n', xaxt='n', yaxt='n')
   legend(...)
 }
 par(mar = c(5, 4, 1.7, 0.2))
 y=as.numeric(format(as.Date(rownames(StockIndex))[60+1:k],"%Y"))
 m=as.numeric(format(as.Date(rownames(StockIndex))[60+1:k],"%m"))
 rownames(P)=y+(m-1)/12
 barplot(t(P)*100,col=darkcols)
 add_legend("topright", legend=colnames( StockIndex ), pch=20, 
            col=darkcols,
            horiz=TRUE, bty='n', cex=1)

or, alternatively

library(FRAPO)
library(fPortfolio)
## Retrieving data and calculating returns
data(StockIndex)
StockReturn <- na.omit(timeSeries(returnseries(StockIndex,
                                               method = "discrete"),
                                  charvec = rownames(StockIndex)))
## Specifying portfolio
pspec <- portfolioSpec()
gmv <-  pspec
end <- time(StockReturn)[60:239]
from <- time(StockReturn)[1:length(end)]
wGMV <- matrix(NA, ncol = ncol(StockReturn), nrow = length(end))
for(i in 1:length(end)){
  series <- window(StockReturn, start = from[i], end = end[i])
  gmvpf <- minvariancePortfolio(data = series, spec = gmv,constraints = "LongOnly")
  wGMV[i, ] <- c(getWeights(gmvpf))
}
P=wGMV[-1,]
rownames(P)=y+(m-1)/12
par(mar = c(5, 4, 1.7, 0.2))
y=as.numeric(format(as.Date(rownames(StockIndex))[60+1:k],"%Y"))
m=as.numeric(format(as.Date(rownames(StockIndex))[60+1:k],"%m"))
barplot(t(P)*100,col=darkcols)
add_legend("topright", legend=colnames( StockIndex ), pch=20, 
           col=darkcols,
           horiz=TRUE, bty='n', cex=1)

Using CVaR risk measure

library(FRAPO)
library(fPortfolio)
data(StockIndex)
pspec <- portfolioSpec()
gmv <-  pspec
cvar <- pspec
setType(cvar) <- "CVaR"
setAlpha(cvar) <- 0.1
setSolver(cvar) <- "solveRglpk.CVAR"
end <- time(StockReturn)[60:239]
from <- time(StockReturn)[1:length(end)]
wGMV <- matrix(NA, ncol = ncol(StockReturn), nrow = length(end))
wCVAR <- wGMV
for(i in 1:length(end)){
  series <- window(StockReturn, start = from[i], end = end[i])
  gmvpf <- minvariancePortfolio(data = series, spec = gmv,constraints = "LongOnly")
  wGMV[i, ] <- c(getWeights(gmvpf))
  cvarpf <- minriskPortfolio(data = series, spec = cvar,constraints = "LongOnly")
  wCVAR[i, ] <- c(getWeights(cvarpf))
}
P=wCVAR[-1,]
rownames(P)=y+(m-1)/12
par(mar = c(5, 4, 1.7, 0.2))
y=as.numeric(format(as.Date(rownames(StockIndex))[60+1:k],"%Y"))
m=as.numeric(format(as.Date(rownames(StockIndex))[60+1:k],"%m"))
barplot(t(P)*100,col=darkcols)
add_legend("topright", legend=colnames( StockIndex ), pch=20, 
           col=darkcols,
           horiz=TRUE, bty='n', cex=1)

Impact on portfolio return

wGMVL1 <- lag(timeSeries(wGMV, charvec = end), k = 1)
colnames(wGMVL1) <- colnames(StockReturn)
wCVARL1 <- lag(timeSeries(wCVAR, charvec = end), k = 1)
colnames(wCVARL1) <- colnames(StockReturn)
GMVRetFac <- 1 + rowSums(wGMVL1 *
                           StockReturn[time(wGMVL1), ]) / 100  
GMVRetFac[1] <- 100
GMVPort <- timeSeries(cumprod(GMVRetFac),
                      charvec = names(GMVRetFac))
CVARRetFac <- 1 + rowSums(wCVARL1 *
                            StockReturn[time(wCVARL1), ]) / 100  
CVARRetFac[1] <- 100
CVARPort <- timeSeries(cumprod(CVARRetFac),
                       charvec = names(CVARRetFac))
ylims <- range(cbind(GMVPort, CVARPort))
plot(GMVPort, ylim = ylims, xlab = "",
     ylab = "Portfolio Value (Index)", col="red")
lines(CVARPort, col = "blue")
legend("topleft",
       legend = c("Minimum-Variance", "Minimum-CVaR"),
       col = c("red", "blue"), lty = 1)

or with constant weights

wGMV0= wGMV
for(i in 2:179) wGMV0[i,]=wGMV0[1,]
wCVAR0=wCVAR
for(i in 2:179) wCVAR0[i,]=wCVAR0[1,]
wGMVL1 <- lag(timeSeries(wGMV0, charvec = end), k = 1)
colnames(wGMVL1) <- colnames(StockReturn)
wCVARL1 <- lag(timeSeries(wCVAR0, charvec = end), k = 1)
colnames(wCVARL1) <- colnames(StockReturn)
GMVRetFac <- 1 + rowSums(wGMVL1 * StockReturn[time(wGMVL1), ]) / 100  
GMVRetFac[1] <- 100
GMVPort <- timeSeries(cumprod(GMVRetFac),charvec = names(GMVRetFac))
CVARRetFac <- 1 + rowSums(wCVARL1 * StockReturn[time(wCVARL1), ]) / 100  
CVARRetFac[1] <- 100
CVARPort <- timeSeries(cumprod(CVARRetFac),charvec = names(CVARRetFac))
ylims <- range(cbind(GMVPort, CVARPort))
plot(GMVPort, ylim = ylims, xlab = "",
     ylab = "Portfolio Value (Index)", col="red")
lines(CVARPort, col = "blue")
legend("topleft",
       legend = c("Minimum-Variance", "Minimum-CVaR"),
       col = c("red", "blue"), lty = 1)

What is we change \(\alpha\) in \(CVaR(\alpha)\) ?

pspec <- portfolioSpec()
gmv <-  pspec
cvar1 <- pspec
cvar2 <- pspec
cvar3 <- pspec

setType(cvar1) <- "CVaR"
setAlpha(cvar1) <- 0.1
setSolver(cvar1) <- "solveRglpk.CVAR"

setType(cvar2) <- "CVaR"
setAlpha(cvar2) <- 0.05
setSolver(cvar2) <- "solveRglpk.CVAR"

setType(cvar3) <- "CVaR"
setAlpha(cvar3) <- 0.01
setSolver(cvar3) <- "solveRglpk.CVAR"

end <- time(StockReturn)[60:239]
from <- time(StockReturn)[1:length(end)]
wGMV <- matrix(NA, ncol = ncol(StockReturn), nrow = length(end))
wCVAR1 <- wCVAR2 <- wCVAR3 <- matrix(NA, ncol = ncol(StockReturn), nrow = length(end))
for(i in 1:length(end)){
  series <- window(StockReturn, start = from[i], end = end[i])
  cvarpf <- minriskPortfolio(data = series, spec = cvar1,constraints = "LongOnly")
  wCVAR1[i, ] <- c(getWeights(cvarpf))
  cvarpf <- minriskPortfolio(data = series, spec = cvar2,constraints = "LongOnly")
  wCVAR2[i, ] <- c(getWeights(cvarpf))
  cvarpf <- minriskPortfolio(data = series, spec = cvar3,constraints = "LongOnly")
  wCVAR3[i, ] <- c(getWeights(cvarpf))
}

wCVARL1 <- lag(timeSeries(wCVAR1, charvec = end), k = 1)
colnames(wCVARL1) <- colnames(StockReturn)
wCVARL2 <- lag(timeSeries(wCVAR2, charvec = end), k = 1)
colnames(wCVARL2) <- colnames(StockReturn)
wCVARL3 <- lag(timeSeries(wCVAR3, charvec = end), k = 1)
colnames(wCVARL3) <- colnames(StockReturn)
CVARRetFac1 <- 1 + rowSums(wCVARL1 * StockReturn[time(wCVARL1), ]) / 100  
CVARRetFac1[1] <- 100
CVARPort1 <- timeSeries(cumprod(CVARRetFac1),charvec = names(CVARRetFac1))
CVARRetFac2 <- 1 + rowSums(wCVARL2 * StockReturn[time(wCVARL2), ]) / 100  
CVARRetFac2[1] <- 100
CVARPort2 <- timeSeries(cumprod(CVARRetFac2),charvec = names(CVARRetFac2))
CVARRetFac3 <- 1 + rowSums(wCVARL3 * StockReturn[time(wCVARL3), ]) / 100  
CVARRetFac3[1] <- 100
CVARPort3 <- timeSeries(cumprod(CVARRetFac3),charvec = names(CVARRetFac3))

ylims <- range(cbind(CVARPort1,CVARPort2,CVARPort3))
plot(CVARPort1, ylim = ylims, xlab = "",
     ylab = "Portfolio Value (Index)", col="blue")
lines(CVARPort2, col = "purple")
lines(CVARPort3, col = "red")
legend("topleft",
       legend = c("Minimum-CVaR 10%","Minimum-CVaR  5%","Minimum-CVaR  1%"),
       col = c("blue","purple","red"), lty = 1)

MultiAsset=StockIndex
pr <- timeSeries(MultiAsset, charvec = rownames(MultiAsset))
data <- returns(pr, methdo = "discrete", percentages = TRUE, trim = TRUE)
NAssets <- ncol(pr)
ANames <- colnames(pr)
Sigma <- cov(data)
mu <- colMeans(data)
hull <- markowitzHull(data, nFrontierPoints = 50)
grid <- feasibleGrid(hull, trace = FALSE)
divers <- bestDiversification(grid, trace = FALSE)
mrc.sd <- function(data, weights){
  Sigma <- cov(data)
  a <- mrc(weights, Sigma)
  sd(a)
}
surf <- riskSurface(divers, FUN = "mrc.sd")
allWeights <- attr(divers, "weights")
idx <- sort(unique(allWeights[, 1]))
dropt <- matrix(0, nrow = length(idx), ncol = 2)
idxRow <- 1:length(idx)
for(j in idx){
  w <- matrix(allWeights[allWeights[, 1] == j, -c(1, 2)], ncol = NAssets)
  divm <- vector()
  length(divm) <- nrow(w)
  for(i in 1:nrow(w)){
    divm[i] <- dr(w[i, ], Sigma)
  }
  divmidx <- which.max(divm)
  wopt <- w[divmidx, ]
  dropt[idxRow[j], ] <- c(crossprod(wopt, mu),
                          sqrt(crossprod(wopt, Sigma) %*% wopt))
}

surfacePlot(surf, type = "filled.contour",
            palette = heat.colors, addHull = TRUE, addGrid = FALSE,
            addAssets = FALSE, xlab = "Target Risk", ylab = "Target Return",
            main = "Convex Hull with Risk Surface:\nStd.Dev. of MRC and MDP-line")
box()
frontier <- portfolioFrontier(data)
MVP <- minvariancePoints(frontier)
TGP <- tangencyPoints(frontier)
sap <- singleAssetPoints(frontier)
wewp <- rep(1/NAssets, NAssets)
mewp <- crossprod(mu, wewp)
sewp <- sqrt(crossprod(wewp, Sigma) %*% wewp)
ERC <- PERC(Sigma)

## Iteration: 0
## pobj: 0
## dobj: 0.793907
## pinf: 1
## dinf: 1
## dgap: 7
## 
## Iteration: 1
## pobj: 0.226732
## dobj: -0.581538
## pinf: 0.169623
## dinf: 0.494202
## dgap: 1.67485
## 
## Iteration: 2
## pobj: 0.373803
## dobj: 0.506046
## pinf: 0.058615
## dinf: 0.158436
## dgap: 0.124398
## 
## Iteration: 3
## pobj: 0.0612196
## dobj: 0.245076
## pinf: 0.0444368
## dinf: 0.0926625
## dgap: 0.0181953
## 
## Iteration: 4
## pobj: -0.255354
## dobj: -0.0780489
## pinf: 0.0387067
## dinf: 0.0536418
## dgap: 0.00461539
## 
## Iteration: 5
## pobj: -0.513268
## dobj: -0.350413
## pinf: 0.0340633
## dinf: 0.0222769
## dgap: 0.000866016
## 
## Iteration: 6
## pobj: -0.515908
## dobj: -0.43796
## pinf: 0.0161835
## dinf: 0.00646059
## dgap: 0.000144857
## 
## Iteration: 7
## pobj: -0.464871
## dobj: -0.449244
## pinf: 0.00323765
## dinf: 0.00098539
## dgap: 1.31597e-05
## 
## Iteration: 8
## pobj: -0.451007
## dobj: -0.449606
## pinf: 0.000290178
## dinf: 8.4799e-05
## dgap: 9.22043e-07
## 
## Iteration: 9
## pobj: -0.449682
## dobj: -0.449609
## pinf: 1.50638e-05
## dinf: 4.38632e-06
## dgap: 4.61978e-08
## 
## Iteration: 10
## pobj: -0.449613
## dobj: -0.449609
## pinf: 7.54433e-07
## dinf: 2.19641e-07
## dgap: 2.31452e-09
## 
## Iteration: 11
## pobj: -0.449609
## dobj: -0.449609
## pinf: 3.77112e-08
## dinf: 1.09789e-08
## dgap: 1.15957e-10
## 
## Optimal solution found.

werc <- Weights(ERC) / 100.0
merc <- crossprod(mu, werc)
serc <- sqrt(crossprod(werc, Sigma) %*% werc)
points(sap, col = "darkgreen", pch = 19, cex = 0.8)
text(sap, ANames, col = "darkred", cex = 0.6, pos = 4)
points(TGP, col = "tan", pch = 19, cex = 2.5)
text(TGP[1], TGP[2], "TGP", col = "purple", cex = 0.5)
points(x = sewp, y = mewp, col = "tan", pch = 19, cex = 2.5)
text(sewp, mewp, "EWP", col="purple", cex = 0.5)
points(x = serc, y = merc, col = "tan", pch = 19, cex = 2.5)
text(serc, merc, "ERC", col="purple", cex = 0.5)
points(MVP, col = "tan", pch = 19, cex = 2.5)
text(MVP[1], MVP[2], "MVP", col = "purple", cex = 0.5)