CoCalc -- EfficientFrontier.R

GitHub Repository: braverock/portfolioanalytics
Path: blob/master/sandbox/paper_analysis/insample/EfficientFrontier.R
¹⁴³³ views
1
setwd("c:/Documents and Settings/Administrator/Desktop/risk budget programs")
2

3
# Equal risk portfolio
4
cAssets = 4;
5
p = priskbudget = 0.95;
6

7
mincriterion = "mES" ; percriskcontribcriterion = "mES";
8

9
# Load programs
10

11
source("R_Allocation/Risk_budget_functions.R"); 
12
library(zoo); library(fGarch); library("PerformanceAnalytics"); library("PortfolioAnalytics")
13

14
clean = TRUE
15
#   Load the data
16

17
firstyear = 1976 ; firstquarter = 1; lastyear = 2010; lastquarter = 2; 
18

19
data = read.table( file= paste("data/","/data.txt",sep="") ,header=T)
20
date = as.Date(data[,1],format="%Y-%m-%d")
21

22

23
monthlyR = indexes = zoo( data[,2:(1+cAssets)] , order.by = date )
24

25
if(clean){
26
indexes = monthlyR = clean.boudt2(monthlyR,alpha=0.05)[[1]]
27
}
28

29
mu = apply(indexes,2,'mean')
30
sigma = cov(indexes)
31
M3 = PerformanceAnalytics:::M3.MM(indexes) 
32
M4 = PerformanceAnalytics:::M4.MM(indexes) 
33

34
# Summary stats individual assets
35

36
head(indexes[,1:2],2); tail(indexes[,1:2],2)
37
apply(indexes[,1:2],2,'mean')*12
38
apply(indexes[,1:2],2,'sd')
39
ES(indexes[,1],method="modified")
40
ES(indexes[,2],method="modified")
41

42
#################################################################################
43
# Make Exhibit 3 Risk budget paper: Efficient frontier plot
44
#################################################################################
45

46
mESfun = function( series ){ return( operMES(  series , alpha = 0.05 , r = 2 ) ) }
47
assetmu = apply( monthlyR , 2 , 'mean' )
48
assetCVaR = apply( monthlyR , 2 , 'mESfun' )
49
minmu = min(assetmu); maxmu = max(assetmu); print(minmu*12); print(maxmu*12);
50

51
# Do the optimization: given a return target minimize the largest percentage CVaR in the portfolio
52
#---------------------------------------------------------------------------------------------------
53

54
# unconstrained solution is Minimum CVaR Concentration portfolio (having the equal risk contribution property):
55
# sol = MinMaxCompCVaRconportfolio(R=monthlyR, Riskupper = Inf ,Returnlower= -Inf )
56
minmu = min( apply(monthlyR,2,'mean' )); 
57
sol = PortfolioOptim(  minriskcriterion = "mES" , MinMaxComp = T, percriskcontribcriterion = "mES" , 
58
                Riskupper = Inf ,  mu = mu, sigma = sigma, M3=M3, M4=M4)
59

60
W               = as.vector( sol[[1]] ) ;   vmu             = as.vector( sol[[2]] )
61
vrisk           = as.vector( sol[[3]] ) ;   mPercrisk       = as.vector( sol[[4]] )  
62
vmaxpercrisk    = max( sol[[4]] )       ;   vmaxriskcontrib = max( sol[[5]] ) 
63

64
for( mutarget in seq(minmu+0.00001,maxmu,0.00001) ){
65
   sol = PortfolioOptim(  minriskcriterion = "mES" , MinMaxComp = T, percriskcontribcriterion = "mES" , 
66
                      Riskupper = Inf ,Returnlower= mutarget,
67
                      mu = mu, sigma = sigma, M3=M3, M4=M4) 
68
   W = rbind( W, as.vector( sol[[1]] ) )
69
   mPercrisk = rbind( mPercrisk , as.vector( sol[[4]] ) )
70
   vmu = c( vmu , as.vector( sol[[2]] ))
71
   vrisk = c( vrisk , as.vector( sol[[3]] ) )
72
   vmaxpercrisk = c( vmaxpercrisk , max( sol[[4]] ) )
73
   vmaxriskcontrib = c( vmaxriskcontrib , max( sol[[5]] ) )
74
}
75

76
# For the highest return targets, a very high penalty parameter is (sometimes) needed
77

78
for( mutarget in c(seq( max(vmu) , maxmu, 0.000001),maxmu) ){
79
   print( c("mutarget equal to",mutarget) )
80
   sol = PortfolioOptim(  minriskcriterion = "mES" , MinMaxComp = T, percriskcontribcriterion = "mES" , 
81
            Riskupper = Inf ,Returnlower= mutarget, penalty = 1e9,
82
            mu = mu, sigma = sigma, M3=M3, M4=M4) 
83
   W = rbind( W, as.vector( sol[[1]] ) )
84
   mPercrisk = rbind( mPercrisk , as.vector( sol[[4]] ) )
85
   vmu = c( vmu , as.vector( sol[[2]] ))
86
   vrisk = c( vrisk , as.vector( sol[[3]] ) )
87
   vmaxpercrisk = c( vmaxpercrisk , max( sol[[4]] ) )
88
   vmaxriskcontrib = c( vmaxriskcontrib , max( sol[[5]] ) )
89
}
90

91

92
if(clean){
93
write.csv( W , file = "EffFrontierMinCVaRConc_weights_clean.csv" )
94
write.csv( mPercrisk , file = "EffFrontierMinCVaRConc_percrisk_clean.csv" )
95
EffFrontier_stats = cbind( vmu , vrisk , vmaxpercrisk , vmaxriskcontrib)
96
write.csv( EffFrontier_stats , file = "EffFrontierMinCVaRConc_stats_clean.csv" )
97
}else{
98
write.csv( W , file = "EffFrontierMinCVaRConc_weights.csv" )
99
write.csv( mPercrisk , file = "EffFrontierMinCVaRConc_percrisk.csv" )
100
EffFrontier_stats = cbind( vmu , vrisk , vmaxpercrisk , vmaxriskcontrib)
101
write.csv( EffFrontier_stats , file = "EffFrontierMinCVaRConc_stats.csv" )
102
}
103

104

105

106
# Do the optimization: given a return target minimize the portfolio CVaR
107
#--------------------------------------------------------------------------
108

109
# unconstrained solution is Minimum CVaR portfolio (having the property that percentage CVaR corresponds to porfolio weights):
110
minmu = min( apply(monthlyR,2,'mean' )); 
111
sol = PortfolioOptim(  minriskcriterion = "mES" , MinMaxComp = F, percriskcontribcriterion = "mES" , 
112
                Riskupper = Inf , mu = mu, sigma = sigma, M3=M3, M4=M4) 
113

114
W               = as.vector( sol[[1]] ) ;   vmu             = as.vector( sol[[2]] )
115
vrisk           = as.vector( sol[[3]] ) ;   mPercrisk       = as.vector( sol[[4]] )  
116
vmaxpercrisk    = max( sol[[4]] )       ;   vmaxriskcontrib = max( sol[[5]] ) 
117

118
for( mutarget in seq(minmu+0.00001,maxmu,0.00001) ){
119
   sol = PortfolioOptim(  minriskcriterion = "mES" , MinMaxComp = F, percriskcontribcriterion = "mES" , 
120
                      Riskupper = Inf ,Returnlower= mutarget,
121
                      mu = mu, sigma = sigma, M3=M3, M4=M4) 
122
   W = rbind( W, as.vector( sol[[1]] ) )
123
   mPercrisk = rbind( mPercrisk , as.vector( sol[[4]] ) )
124
   vmu = c( vmu , as.vector( sol[[2]] ))
125
   vrisk = c( vrisk , as.vector( sol[[3]] ) )
126
   vmaxpercrisk = c( vmaxpercrisk , max( sol[[4]] ) )
127
   vmaxriskcontrib = c( vmaxriskcontrib , max( sol[[5]] ) )
128
}
129

130
# For the highest return targets, a very high penalty parameter is (sometimes) needed
131

132
for( mutarget in c(seq( max(vmu) , maxmu, 0.000001),maxmu) ){
133
   print( c("mutarget equal to",mutarget) )
134
   sol = PortfolioOptim(  minriskcriterion = "mES" , MinMaxComp = F, percriskcontribcriterion = "mES" , 
135
            Riskupper = Inf ,Returnlower= mutarget, penalty = 1e9,
136
            mu = mu, sigma = sigma, M3=M3, M4=M4) 
137
   W = rbind( W, as.vector( sol[[1]] ) )
138
   mPercrisk = rbind( mPercrisk , as.vector( sol[[4]] ) )
139
   vmu = c( vmu , as.vector( sol[[2]] ))
140
   vrisk = c( vrisk , as.vector( sol[[3]] ) )
141
   vmaxpercrisk = c( vmaxpercrisk , max( sol[[4]] ) )
142
   vmaxriskcontrib = c( vmaxriskcontrib , max( sol[[5]] ) )
143
}
144

145
if(clean){
146
   write.csv( W , file = "EffFrontierMinCVaR_weights_clean.csv" )
147
   write.csv( mPercrisk , file = "EffFrontierMinCVaR_percrisk_clean.csv" )
148
   EffFrontier_stats = cbind( vmu , vrisk , vmaxpercrisk , vmaxriskcontrib)
149
   write.csv( EffFrontier_stats , file = "EffFrontierMinCVaR_stats_clean.csv" )
150
}else{
151
   write.csv( W , file = "EffFrontierMinCVaR_weights.csv" )
152
   write.csv( mPercrisk , file = "EffFrontierMinCVaR_percrisk.csv" )
153
   EffFrontier_stats = cbind( vmu , vrisk , vmaxpercrisk , vmaxriskcontrib)
154
   write.csv( EffFrontier_stats , file = "EffFrontierMinCVaR_stats.csv" )
155
}
156

157

158

159

160
# Do the optimization: given a return target minimize the portfolio standard deviation (use quadprog)
161
#------------------------------------------------------------------------------------------------------
162

163
library(quadprog)
164
N = 4 ;  maxweight = 1; minweight = 0; sumweights = 1; 
165
Amat = rbind( rep(1,N) ,  diag(x =1,nrow=N,ncol=N) ,  diag(x =-1,nrow=N,ncol=N) , as.numeric(mu) );
166
bvec =  c( sumweights , rep(minweight,N), rep(-maxweight,N)  ) 
167
dvec = matrix( rep(0,N) , ncol=1 )
168

169
mutarget = -10000; 
170
# min(-d^T b + 1/2 b^T D b) with the constraints A^T b >= b_0.
171
optw = solve.QP( Dmat = sigma   , dvec = dvec , Amat = t(Amat) , 
172
   bvec = matrix( c(bvec,mutarget) , ncol = 1) , meq =1  )$solution
173
sol = ES(weights=optw, portfolio_method="component", mu = mu, sigma = sigma, m3=M3, m4=M4,invert=FALSE)
174
W               = optw                             ;   vmu             = sum(optw*mu)
175
vrisk           = sol$MES                          ;   mPercrisk       = as.vector( sol$pct_contrib_MES )  
176
vmaxpercrisk    = max( sol$pct_contrib_MES )       ;   vmaxriskcontrib = max( sol$contribution ) 
177

178
# unconstrained solution is Minimum CVaR portfolio (having the property that percentage CVaR corresponds to porfolio weights):
179
minmu = min( apply(monthlyR,2,'mean' )); 
180

181

182
for( mutarget in seq(minmu+0.00001,maxmu,0.00001) ){
183
   optw = solve.QP( Dmat = sigma   , dvec = dvec , Amat = t(Amat) , 
184
          bvec = matrix( c(bvec,mutarget) , ncol = 1) , meq =1  )$solution
185
   sol = ES(weights=optw, portfolio_method="component", mu = mu, sigma = sigma, m3=M3, m4=M4,invert=FALSE)
186
   W = rbind( W, optw )
187
   mPercrisk = rbind( mPercrisk , as.vector( sol$pct_contrib_MES ) )
188
   vmu = c( vmu , sum(optw*mu) )
189
   vrisk = c( vrisk , sol$MES )
190
   vmaxpercrisk = c( vmaxpercrisk , max( sol$pct_contrib_MES ) )
191
   vmaxriskcontrib = c( vmaxriskcontrib , max( sol$contribution ) )
192
}
193

194
if(clean){
195
   write.csv( W , file = "EffFrontierMinVar_weights_clean.csv" )
196
   write.csv( mPercrisk , file = "EffFrontierMinVar_percrisk_clean.csv" )
197
   EffFrontier_stats = cbind( vmu , vrisk , vmaxpercrisk , vmaxriskcontrib)
198
   write.csv( EffFrontier_stats , file = "EffFrontierMinVar_stats_clean.csv" )
199
}else{
200
   write.csv( W , file = "EffFrontierMinVar_weights.csv" )
201
   write.csv( mPercrisk , file = "EffFrontierMinVar_percrisk.csv" )
202
   EffFrontier_stats = cbind( vmu , vrisk , vmaxpercrisk , vmaxriskcontrib)
203
   write.csv( EffFrontier_stats , file = "EffFrontierMinVar_stats.csv" )
204
}
205

206

207

208

209
# Plotting 
210
#----------------------------------------------------------------------
211

212
# Layout 3 
213
source("R_interpretation/chart.StackedBar.R"); 
214

215
mESfun = function( series ){ return( operMES(  series , alpha = 0.05 , r = 2 ) ) }
216
assetmu = apply( monthlyR , 2 , 'mean' )
217
assetCVaR = apply( monthlyR , 2 , 'mESfun' )
218
minmu = min(assetmu); maxmu = max(assetmu); print(minmu*12); print(maxmu*12);
219

220

221

222
####################################
223
# Min CVaR Concentration
224
####################################
225
labelnames = c( "US bond" , "S&P 500" , "EAFE" , "GSCI" )
226

227
clean = TRUE
228
if(clean){
229
 # Portfolio weights
230
 W_MCC                 = read.csv( file = "EffFrontierMinCVaRConc_weights_clean.csv" )[,2:(1+cAssets)]
231
 W_minCVaR             = read.csv( file = "EffFrontierMinCVaR_weights_clean.csv" )[,2:(1+cAssets)]
232
 W_minVar              = read.csv( file = "EffFrontierMinVar_weights_clean.csv" )[,2:(1+cAssets)]
233
 # Percentage CVaR contributions
234
 PercCVaR_MCC         = read.csv(  file = "EffFrontierMinCVaRConc_percrisk_clean.csv" )[,2:(1+cAssets)]
235
 PercCVaR_minCVaR     = read.csv(  file = "EffFrontierMinCVaR_percrisk_clean.csv" )[,2:(1+cAssets)]
236
 PercCVaR_minVar      = read.csv(  file = "EffFrontierMinVar_percrisk_clean.csv" )[,2:(1+cAssets)]
237
 # Summary stats
238
 EffFrontier_stats_MCC = read.csv(file = "EffFrontierMinCVaRConc_stats_clean.csv")[,2:5]
239
 EffFrontier_stats_minVar = read.csv(file = "EffFrontierMinVar_stats_clean.csv")[,2:5]
240
 EffFrontier_stats_minCVaR = read.csv(file = "EffFrontierMinCVar_stats_clean.csv")[,2:5]
241
}else{
242
 # Portfolio weights
243
 W_MCC                 = read.csv(file = "EffFrontierMinCVaRConc_weights.csv")[,2:(1+cAssets)]
244
 W_minCVaR             = read.csv(file = "EffFrontierMinCVaR_weights.csv")[,2:(1+cAssets)]
245
 W_minVar              = read.csv(file = "EffFrontierMinVar_weights.csv")[,2:(1+cAssets)]
246
 # Percentage CVaR contributions
247
 mPercrisk_MCC         = read.csv( file = "EffFrontierMinCVaRConc_percrisk.csv" )[,2:(1+cAssets)]
248
 mPercrisk_minCVaR     = read.csv( file = "EffFrontierMinCVaR_percrisk.csv" )[,2:(1+cAssets)]
249
 mPercrisk_minVar      = read.csv(file = "EffFrontierMinVar_percrisk.csv")[,2:(1+cAssets)]
250
 # Summary stats
251
 EffFrontier_stats_MCC       = read.csv(file = "EffFrontierMinCVaRConc_stats.csv")[,2:5]
252
 EffFrontier_stats_minVar = read.csv(file = "EffFrontierMinVar_stats.csv")[,2:5]
253
 EffFrontier_stats_minCVaR   = read.csv(file = "EffFrontierMinCVar_stats.csv")[,2:5]
254
}
255

256
vmu_MCC             = EffFrontier_stats_MCC[,1] ;  vmu_minCVaR             = EffFrontier_stats_minCVaR[,1] ;
257
vrisk_MCC           = EffFrontier_stats_MCC[,2] ;  vrisk_minCVaR           = EffFrontier_stats_minCVaR[,2] ;
258
vmaxpercrisk_MCC    = EffFrontier_stats_MCC[,3] ;  vmaxpercrisk_minCVaR    = EffFrontier_stats_minCVaR[,3] ;
259
vriskconc_MCC       = EffFrontier_stats_MCC[,4] ;  vriskconc_minCVaR       = EffFrontier_stats_minCVaR[,4] ;
260

261
vmu_minVar        = EffFrontier_stats_minVar[,1] ;
262
vrisk_minVar      = EffFrontier_stats_minVar[,2] ;
263
vriskconc_minVar  = EffFrontier_stats_minVar[,4] ;
264

265
# The MCC and Min CVaR were obtained with DEoptim. The solutions are not always optimal. Correct for this: 
266

267
order_MCC = sort(vmu_MCC,index.return=T)$ix      ; order_minCVaR = sort(vmu_minCVaR,index.return=T)$ix 
268
vmu_MCC = vmu_MCC[order_MCC]                     ; vmu_minCVaR = vmu_minCVaR[order_minCVaR] 
269
vrisk_MCC = vrisk_MCC[order_MCC]                 ; vrisk_minCVaR = vrisk_minCVaR[order_minCVaR]      
270
vmaxpercrisk_MCC = vmaxpercrisk_MCC[order_MCC]   ; vmaxpercrisk_minCVaR = vmaxpercrisk_MCC[order_minCVaR] 
271
vriskconc_MCC = vriskconc_MCC[order_MCC]         ; vriskconc_minCVaR = vriskconc_minCVaR[order_minCVaR] 
272
W_MCC  = W_MCC[order_MCC,]                       ; W_minCVaR  = W_minCVaR[order_minCVaR,]     
273
PercCVaR_MCC = PercCVaR_MCC[order_MCC,]          ; PercCVaR_minCVaR = PercCVaR_minCVaR[order_minCVaR,]   
274

275
# eliminate duplicates in mu
276
# Determines which elements of a vector or data frame are duplicates of elements with smaller subscripts, 
277
# and returns a logical vector indicating which elements (rows) are duplicates.
278
a_MCC = duplicated(vmu_MCC)                         ; a_minCVaR = duplicated(vmu_minCVaR)  
279
vmu_MCC = vmu_MCC[!a_MCC]                           ; vmu_minCVaR = vmu_minCVaR[!a_minCVaR] 
280
vrisk_MCC = vrisk_MCC[!a_MCC]                       ; vrisk_minCVaR = vrisk_minCVaR[!a_minCVaR]
281
vmaxpercrisk_MCC = vmaxpercrisk_MCC[!a_MCC]         ; vmaxpercrisk_minCVaR = vmaxpercrisk_minCVaR[!a_minCVaR] 
282
vriskconc_MCC = vriskconc_MCC[!a_MCC]               ; vriskconc_minCVaR = vriskconc_minCVaR[!a_minCVaR]  
283
W_MCC = W_MCC[!a_MCC,]                              ; W_minCVaR = W_minCVaR[!a_minCVaR,]    
284
PercCVaR_MCC = PercCVaR_MCC[!a_MCC,]                ; PercCVaR_minCVaR = PercCVaR_minCVaR[!a_minCVaR,] 
285
 
286
# eliminate (efficiency: lower mu, higher CVaR alloc
287

288
sel_MCC = sel_minCVaR = c(); nmu_MCC = length(vmu_MCC); nmu_minCVaR = length(vmu_minCVaR); 
289
for( i in 1:(nmu_MCC-1) ){
290
    if( any(vriskconc_MCC[(i+1):nmu_MCC]<=vriskconc_MCC[i]) ){ }else{sel_MCC = c(sel_MCC, i) }
291
}
292
for( i in 1:(nmu_minCVaR-1) ){
293
    if( any(vrisk_minCVaR[(i+1):nmu_minCVaR]<=vrisk_minCVaR[i]) ){ }else{sel_minCVaR = c(sel_minCVaR, i) }
294
}
295
vmu_MCC = vmu_MCC[sel_MCC]                        ; vmu_minCVaR = vmu_minCVaR[sel_minCVaR]  
296
vrisk_MCC = vrisk_MCC[sel_MCC]                    ; vrisk_minCVaR = vrisk_minCVaR[sel_minCVaR] 
297
vmaxpercrisk_MCC = vmaxpercrisk_MCC[sel_MCC]      ; vmaxpercrisk_minCVaR = vmaxpercrisk_minCVaR[sel_minCVaR]
298
vriskconc_MCC = vriskconc_MCC[sel_MCC]            ; vriskconc_minCVaR = vriskconc_minCVaR[sel_minCVaR] 
299
W_MCC = W_MCC[sel_MCC,]                           ; W_minCVaR = W_minCVaR[sel_minCVaR,]
300
PercCVaR_MCC = PercCVaR_MCC[sel_MCC,]             ; PercCVaR_minCVaR = PercCVaR_minCVaR[sel_minCVaR,] 
301

302
wEW <- rep(1/4,4)
303
muEW       = mean(assetmu*12)
304
outES       = ES(weights=wEW, portfolio_method="component", mu = mu, sigma = sigma, m3=M3, m4=M4,invert=FALSE) 
305
riskEW     = outES$MES
306
riskconcEW = max(outES$contribution)
307

308
postscript('frontier_fourassets.eps') 
309

310
layout( matrix(  c(1,2,3,1,2,3),  ncol = 2 ) , height= c(5,5,1), width=c(1,1) )
311
ylim = c( min(assetmu) , max(assetmu) )
312
par( mar = c(4.5,5,2,2), las=1 ,cex=0.9 , cex.axis=0.9, cex.main=0.95)
313
plot( vrisk_MCC , vmu_MCC*12 , type="l", lty = 1, 
314
          main = "" , 
315
          ylab="Annualized mean return" , xlab="Monthly 95% Portfolio CVaR" , 
316
          lwd=2, xlim = c(0.02,0.13) , ylim = (ylim*12+c(0,0.01)) , col="black" )
317
lines( vrisk_minCVaR , vmu_minCVaR*12,  lty = 1, lwd=2, col="darkgray" ) 
318
lines( vrisk_minVar , vmu_minVar*12,  lty = 3, lwd=2, col="black" )
319

320
c = 1; text(y=assetmu[c]*12,x= assetCVaR[c]  , label=labelnames[c]  , cex = 0.7, offset = 0.2, pos = 1); points(y=assetmu[c]*12,x= assetCVaR[c] )  
321
for( c in 2:cAssets ){ 
322
   text(y=assetmu[c]*12,x= assetCVaR[c]  , label=labelnames[c]  , cex = 0.7, offset = 0.2, pos = 3) 
323
   points(y=assetmu[c]*12,x= assetCVaR[c] ) 
324
};
325
# Plot also EW portfolio
326

327
text(y=muEW,x=riskEW, label="EW"  , cex = 0.7, offset = 0.2, pos = 3)
328
points(y=muEW,x= riskEW ) 
329

330
par( mar = c(4.5,5,2,2), las=1 ,cex=0.9 , cex.axis=0.9, cex.main=0.95)
331
plot( vriskconc_MCC/vrisk_MCC , vmu_MCC*12 , type="l", lty = 1, 
332
          main = "" , 
333
          ylab="Annualized  mean return" , xlab="Largest Perc. CVaR contribution" , 
334
          lwd=2, xlim = c(0.2,1.1), ylim = c(ylim*12+c(0,0.01)) , col="black" )
335
lines( vriskconc_minCVaR/vrisk_minCVaR , vmu_minCVaR*12,  lty = 1, lwd=2, col="darkgray" ) 
336
lines( vriskconc_minVar/vrisk_minVar , vmu_minVar*12,  lty = 3, lwd=2, col="black" )
337

338
text(y=muEW,x=riskconcEW/riskEW, label="EW"  , cex = 0.7, offset = 0.2, pos = 1)
339
points(y=muEW,x= riskconcEW/riskEW ) 
340
for( c in 1:cAssets ){ 
341
   text(y=assetmu[c]*12,x= 1  , label=labelnames[c]  , cex = 0.7, offset = 0.2, pos = 4) 
342
   points(y=assetmu[c]*12,x= 1 ) 
343
};
344

345
par( mar = c(0,1,1,0) )
346
plot.new()
347
legend("center",legend=c("Mean/StdDev" , "Mean/CVaR","Mean/CVaR concentration"  ),
348
   lty=c("dashed","solid","solid"), cex=0.8,ncol=3,lwd=2,col=c("black","darkgray","black"))
349
dev.off()
350

351

352
# cbind( vmu_minVar*12 , vriskconc_minVar/vrisk_minVar ) : 0.0816 is turning point
353
# cbind( vmu_minCVaR*12 , vriskconc_minCVaR/vrisk_minCVaR ) : 0.0816 is turning point
354

355

356
# Make the weight/CVaR allocation plots
357
#---------------------------------------
358

359
Wsel_MCC = PercCVaRsel_MCC  = Wsel_minCVaR = PercCVaRsel_minCVaR  = Wsel_minVar = PercCVaRsel_minVar  = c(); 
360
vmusel_MCC = vrisksel_MCC = vriskconcsel_MCC = vmusel_minCVaR = vrisksel_minCVaR = vriskconcsel_minCVaR =  vmusel_minVar = vrisksel_minVar = vriskconcsel_minVar =   c();
361
lm = minmu = min( c(vmu_MCC,vmu_minCVaR) ) ; maxmu = max( c(vmu_MCC,vmu_minCVaR) );
362
ylim = c( min(assetmu) - 0.0003 , max(assetmu) + 0.0003 )
363
xlim = c( 0 , max(assetCVaR) + 0.01 )
364

365
binning = T # for the weight plots, binned data such that no visual misinterpretation is possible
366
step = quantile( diff(vmu_MCC) , 0.1 ) 
367
if( binning ){
368
   for( rm in seq( minmu+step , maxmu , step )  ){
369
      selection    =  c(vmu_MCC >= lm & vmu_MCC < rm) ;
370
      if( any(selection) ){
371
         selection    = c(1:length(selection))[selection]
372
         selone       = sort(vriskconc_MCC[ selection ],index.return=T)$ix[1]
373
         selone       = selection[selone]
374
         vmusel_MCC       = c( vmusel_MCC       ,  mean(vmu_MCC[selone  ] ));
375
         Wsel_MCC         = rbind( Wsel_MCC     ,  apply( W_MCC[selone,] ,2,'mean' ))
376
         PercCVaRsel_MCC  = rbind( PercCVaRsel_MCC  ,  apply( PercCVaR_MCC[selone,] ,2,'mean' )) 
377
         vrisksel_MCC     = c( vrisksel_MCC     ,  mean(vrisk_MCC[ selone ]) );
378
         vriskconcsel_MCC = c( vriskconcsel_MCC ,  mean(vriskconc_MCC[ selone ])  )
379
      }else{
380
         vmusel_MCC       = c( vmusel_MCC       ,  NA );
381
         Wsel_MCC         = rbind( Wsel_MCC     ,  rep(NA,cAssets)  )
382
         PercCVaRsel_MCC  = rbind( PercCVaRsel_MCC  ,  rep(NA,cAssets) ) 
383
         vrisksel_MCC     = c( vrisksel_MCC     ,  NA );
384
         vriskconcsel_MCC = c( vriskconcsel_MCC ,  NA  )
385
      }
386
      selection    =  c(vmu_minCVaR >= lm & vmu_minCVaR < rm) ;
387
      if( any(selection) ){
388
         selection    = c(1:length(selection))[selection]
389
         selone       = sort(vrisk_minCVaR[ selection ],index.return=T)$ix[1]
390
         selone       = selection[selone]
391
         vmusel_minCVaR       = c( vmusel_minCVaR       ,  mean(vmu_minCVaR[selone  ] ));
392
         Wsel_minCVaR         = rbind( Wsel_minCVaR     ,  apply( W_minCVaR[selone,] ,2,'mean' ))
393
         PercCVaRsel_minCVaR  = rbind( PercCVaRsel_minCVaR  ,  apply( PercCVaR_minCVaR[selone,] ,2,'mean' )) 
394
         vrisksel_minCVaR     = c( vrisksel_minCVaR     ,  mean(vrisk_minCVaR[ selone ]) );
395
         vriskconcsel_minCVaR = c( vriskconcsel_minCVaR ,  mean(vriskconc_minCVaR[ selone ])  )
396
      }else{
397
         vmusel_minCVaR       = c( vmusel_minCVaR       ,  NA );
398
         Wsel_minCVaR         = rbind( Wsel_minCVaR     ,  rep(NA,cAssets) )
399
         PercCVaRsel_minCVaR  = rbind( PercCVaRsel_minCVaR  ,  rep(NA,cAssets) ) 
400
         vrisksel_minCVaR     = c( vrisksel_minCVaR     ,  NA );
401
         vriskconcsel_minCVaR = c( vriskconcsel_minCVaR ,  NA )
402
      }
403
      selection    =  c(vmu_minVar >= lm & vmu_minVar < rm) ;
404
      if( any(selection) ){
405
         selection    = c(1:length(selection))[selection]
406
         selone       = sort(vrisk_minVar[ selection ],index.return=T)$ix[1]
407
         selone       = selection[selone]
408
         vmusel_minVar       = c( vmusel_minVar       ,  mean(vmu_minVar[selone  ] ));
409
         Wsel_minVar         = rbind( Wsel_minVar     ,  apply( W_minVar[selone,] ,2,'mean' ))
410
         PercCVaRsel_minVar  = rbind( PercCVaRsel_minVar  ,  apply( PercCVaR_minVar[selone,] ,2,'mean' )) 
411
         vrisksel_minVar     = c( vrisksel_minVar     ,  mean(vrisk_minVar[ selone ]) );
412
         vriskconcsel_minVar = c( vriskconcsel_minVar ,  mean(vriskconc_minVar[ selone ])  )
413
      }else{
414
         vmusel_minVar       = c( vmusel_minVar       ,  NA );
415
         Wsel_minVar         = rbind( Wsel_minVar     ,  rep(NA,cAssets) )
416
         PercCVaRsel_minVar  = rbind( PercCVaRsel_minVar  ,  rep(NA,cAssets) ) 
417
         vrisksel_minVar     = c( vrisksel_minVar     ,  NA );
418
         vriskconcsel_minVar = c( vriskconcsel_minVar ,  NA  )
419
      }
420

421
      lm = rm; 
422
   }
423
}else{
424
   vmusel_MCC = vmu_MCC              ;   vmusel_minCVaR = vmu_minCVaR 
425
   Wsel_MCC   = W_MCC                ;   Wsel_minCVaR   = W_minCVaR                ;
426
   vrisksel_MCC = vrisk_MCC          ;   vrisksel_minCVaR = vrisk_minCVaR  
427
   vriskconcsel_MCC = vriskconc_MCC  ;   vriskconcsel_minCVar = vriskconc_minCVaR
428
   PercCVaRsel_MCC =  PercCVaR_MCC   ;   PercCVaRsel_minCVar =  PercCVaR_minCVar       ; 
429
}
430

431
library(zoo)
432
cAssets = ncol(monthlyR)
433
 colorset = gray( seq(0,(cAssets-1),1)/cAssets ) ; 
434
colnames( Wsel_MCC ) = colnames( Wsel_minCVaR ) = labelnames
435
rownames( Wsel_MCC ) = rownames( PercCVaRsel_MCC ) = round( interpNA(vmusel_MCC)*12,4); 
436
rownames( Wsel_minCVaR ) = rownames( PercCVaRsel_minCVaR ) = round( interpNA(vmusel_minCVaR)*12,4) ;
437

438
rownames( Wsel_minVar ) = rownames( PercCVaRsel_minVar ) = round( interpNA(vmusel_minVar)*12,4) ; 
439
colnames( Wsel_minVar ) = labelnames
440

441

442

443
w.names = c( "US bond" , "S&P 500", "EAFE"  , "GSCI" )
444
namelabels = c("Mean/StdDev" , "Mean/CVaR","Mean/CVaR concentration"  )
445
 l = 2
446
mar1 =c(2,l,2,1.1)
447
mar2 =c(0,l,2,1)
448
mar3 = c(3,l+1,3,0.1)
449
mar4 = c(2,l+1,2,0.1)
450

451
source("R_interpretation/chart.StackedBar.R"); 
452

453
# Stacked weights plot: 
454
postscript('stackedweightsriskcont_efficientfrontier_withNA.eps') 
455
   layout( matrix(  c(1,2,3,4,5,6,7,4),  ncol = 2 ) , height= c(1.5,1.5,1.5,0.7), width=1)
456

457
   par(mar=mar3 , cex.main=1)
458
   chart.StackedBar2(Wsel_minVar    ,col=colorset,space=0,  main = namelabels[1], ylab="Weight allocation", las=1, l=3.9, r=0,  cex.axis=1, cex.lab=1,  cex.main=1, axisnames=T,legend.loc = NULL,ylim=c(0,1),border = F )
459
   chart.StackedBar2(Wsel_minCVaR,col=colorset,space=0,  main = namelabels[2], ylab="Weight allocation", las=1, l=3.9, r=0,   cex.axis=1, cex.lab=1,  cex.main=1, axisnames=T,legend.loc = NULL,ylim=c(0,1),border = F )
460
   chart.StackedBar2(Wsel_MCC    ,col=colorset,space=0,  main = namelabels[3], ylab="Weight allocation", las=1, l=3.9, r=0, cex.axis=1, cex.lab=1,  cex.main=1, axisnames=T, legend.loc = NULL,ylim=c(0,1),border = F )
461

462
   par(mar=mar1 , cex.main=1)
463
   plot.new()
464
   legend("center",legend=w.names,fill=colorset,ncol=4)
465

466
   par(mar=mar3 , cex.main=1)
467
   chart.StackedBar2(PercCVaRsel_minVar,col=colorset,space=0,  main = namelabels[1], ylab="CVaR allocation", las=1, l=3.9, r=0,  cex.axis=1, cex.lab=1,  cex.main=1, axisnames=T,legend.loc = NULL,ylim=c(0,1),border = F )
468
   chart.StackedBar2(PercCVaRsel_minCVaR,col=colorset,space=0,  main = namelabels[2], ylab="CVaR allocation", las=1, l=3.9, r=0,   cex.axis=1, cex.lab=1,  cex.main=1, axisnames=T,legend.loc = NULL,ylim=c(0,1),border = F )
469
   chart.StackedBar2(PercCVaRsel_MCC,col=colorset,space=0,  main = namelabels[3], ylab="CVaR allocation", las=1, l=3.9, r=0,  cex.axis=1, cex.lab=1,  cex.main=1, axisnames=T,legend.loc = NULL,ylim=c(0,1),border = F )
470

471
dev.off()
472

473

474
postscript('stackedweightsriskcont_efficientfrontier.eps') 
475
   layout( matrix(  c(1,2,3,4,5,6,7,4),  ncol = 2 ) , height= c(1.5,1.5,1.5,0.7), width=1)
476

477
   par(mar=mar3 , cex.main=1)
478
   chart.StackedBar2(interpNA(Wsel_minVar)    ,col=colorset,space=0,  main = namelabels[1], ylab="Weight allocation", las=1, l=3.9, r=0,  cex.axis=1, cex.lab=1,  cex.main=1, axisnames=T,legend.loc = NULL,ylim=c(0,1),border = F )
479
   chart.StackedBar2(interpNA(Wsel_minCVaR),col=colorset,space=0,  main = namelabels[2], ylab="Weight allocation", las=1, l=3.9, r=0,   cex.axis=1, cex.lab=1,  cex.main=1, axisnames=T,legend.loc = NULL,ylim=c(0,1),border = F )
480
   chart.StackedBar2(interpNA(Wsel_MCC)    ,col=colorset,space=0,  main = namelabels[3], ylab="Weight allocation", las=1, l=3.9, r=0, cex.axis=1, cex.lab=1,  cex.main=1, axisnames=T, legend.loc = NULL,ylim=c(0,1),border = F )
481

482
   par(mar=mar1 , cex.main=1)
483
   plot.new()
484
   legend("center",legend=w.names,fill=colorset,ncol=4)
485

486
   par(mar=mar3 , cex.main=1)
487
   chart.StackedBar2(interpNA(PercCVaRsel_minVar),col=colorset,space=0,  main = namelabels[1], ylab="CVaR allocation", las=1, l=3.9, r=0,  cex.axis=1, cex.lab=1,  cex.main=1, axisnames=T,legend.loc = NULL,ylim=c(0,1),border = F )
488
   chart.StackedBar2(interpNA(PercCVaRsel_minCVaR),col=colorset,space=0,  main = namelabels[2], ylab="CVaR allocation", las=1, l=3.9, r=0,   cex.axis=1, cex.lab=1,  cex.main=1, axisnames=T,legend.loc = NULL,ylim=c(0,1),border = F )
489
   chart.StackedBar2(interpNA(PercCVaRsel_MCC),col=colorset,space=0,  main = namelabels[3], ylab="CVaR allocation", las=1, l=3.9, r=0,  cex.axis=1, cex.lab=1,  cex.main=1, axisnames=T,legend.loc = NULL,ylim=c(0,1),border = F )
490

491
dev.off()
492

493

494

495

496
#----------------------------------------------------------------------------------------------------
497
# DISCUSSION RESULTS
498
# On Raw Data: Discontinuity in the Mean-CVaR frontier of Mean CVaR concentration efficient portfolios: 
499
# Two types of portfolios:
500
# Lower return targets: portfolio where the component risk of the bond dominates
501
# Higher return targets: portfolio where the component risk of the US equity dominates
502
# Shift when: 
503
# > vmu[13:14]
504
# [1] 0.006796657 0.006806883
505
# > vriskconc[13:14]
506
# [1] 0.0082 0.0254
507
#  > vrisk[13:14]
508
# [1] 0.01544720 0.06036618
509

510
# On Cleaned Data: Three segments
511

512
# > assetmu
513
#      Bond      SP500       EAFE     SPGSCI 
514
# 0.07545960 0.10249338 0.08677129 0.05413622 
515
# > assetCVaR
516
#      Bond      SP500       EAFE     SPGSCI 
517
# 0.02456233 0.09974267 0.10868066 0.12775266 
518

519

520
# Unconstrained: Equal Risk Contribution Portfolio:
521
# > W[1,]
522
# 0.6431897 0.1337981 0.1044994 0.1185127
523
#> mPercrisk[1,]
524
# 0.2516 0.2517 0.2504 0.2463
525
# > vmu[1]*12
526
# [1] 0.07773165
527
#> vrisk[1]
528
#[1] 0.03467514
529
#> vriskconc[1]
530
#[1] 0.0087
531
# Segment 1: increasing allocation to bonds and decreasing allocation to commodities
532
# Portfolio risk concentration increases, portfolio risk decreases
533
# > W[25:27,]
534
# 0.7052788 0.1584357 0.1321813 0.0041042550
535
# 0.7085301 0.1597522 0.1313759 0.0003417814
536
# 0.7025028 0.1628656 0.1344935 0.0001380022
537
#> vmu[26]*12
538
#[1] 0.0812571
539
#> vrisk[26]
540
#[1] 0.03352778
541
#> vriskconc[26]
542
#[1] 0.0112
543
# Segment 2: increasing allocation to both US equity and EAFE
544
# Portfolio risk concentration and risk increase together
545
# > W[138:140,]
546
# 0.0078323664 0.5106416 0.4814090 1.169711e-04
547
# 0.0005986683 0.5136429 0.4856355 1.228907e-04
548
# 0.0002258733 0.5201106 0.4796231 4.051641e-05
549
# > mPercrisk[138:140,]
550
#  0e+00 0.4999 0.5000  0
551
#> vmu[139]*12
552
#[1] 0.09483605
553
# > vrisk[139]
554
# [1] 0.09808784
555
# > vriskconc[139]
556
# [1] 0.049
557
# Segment 3: increasing alllocation to US equity and decreasing to EAFE
558
# Portfolio risk concentration and risk increase together
559
# > tail(W,1)
560
# 5.45402e-05 0.9998762 4.340072e-05 2.585630e-05
561

562

563

564
Product

Resources

Company
