stackedweightriskcontributionplot.R

#################################################################################
# Create stacked weights and risk contributions plot
#################################################################################

# ! Set your working directory (folder containing the subfolders R_allocation, R_interpretation, data, weights, etc)

setwd("c:/Documents and Settings/Administrator/Desktop/risk budget programs")


# Options:
################

# specify the number of years used for the estimation
estyears = 8;

# Load programs

source("R_interpretation/chart.StackedBar.R"); 
library(zoo);  library(PerformanceAnalytics)

# number of risky assets
firstyear = 1976 ; firstquarter = 1; lastyear = 2010; lastquarter = 2; 
cAssets = 4 
 
# "MinRisk_ReturnTarget"  "EqualRisk"
# "MinRiskConc_ReturnTarget"

names = c(  "EqualWeight" , "MinRisk" , "MinRisk_PositionLimit" , "MinRisk_RiskLimit" ,
            "MinRiskConc" , "MinRiskConc_PositionLimit", "EqualRisk" ,
            "MinRisk_ReturnTarget", "MinRiskConc_ReturnTarget" )

namelabels = c( "Equal-Weight" , "Min CVaR" , "Min CVaR + 40% Position Limit" , "Min CVaR + 40% CVaR Alloc Limit" ,  
       "Min CVaR Concentration" , "Min CVaR Concentration + 40% Position Limit", "Min CVaR + ERC constraint" , "Min CVaR + Return Target" , "Min CVaR Conc + Return Target" )

# frequency of rebalancing: yearly of quarterly
frequency = "quarterly"
# Load portfolio weights:
weightsS1 = read.csv( file = paste("weights/", names[1], ".csv" , sep="")  );
weightsS2 = read.csv( file = paste("weights/", names[2], ".csv" , sep="")  );
weightsS3 = read.csv( file = paste("weights/", names[3], ".csv" , sep="")  );
weightsS4 = read.csv( file = paste("weights/", names[4], ".csv" , sep="")  );
weightsS5 = read.csv( file = paste("weights/", names[5], ".csv" , sep="")  );
weightsS6 = read.csv( file = paste("weights/", names[6], ".csv" , sep="")  );
weightsS7 = read.csv( file = paste("weights/", names[7], ".csv" , sep="")  );


# Load percentage risk contributions:
riskcontS1 = read.csv( file = paste("riskcont/", names[1], ".csv" , sep="")  );
riskcontS2 = read.csv( file = paste("riskcont/", names[2], ".csv" , sep="")  );
riskcontS3 = read.csv( file = paste("riskcont/", names[3], ".csv" , sep="")  );
riskcontS4 = read.csv( file = paste("riskcont/", names[4], ".csv" , sep="")  );
riskcontS5 = read.csv( file = paste("riskcont/", names[5], ".csv" , sep="")  );
riskcontS6 = read.csv( file = paste("riskcont/", names[6], ".csv" , sep="")  );
riskcontS7 = read.csv( file = paste("riskcont/", names[7], ".csv" , sep="")  );

if(frequency=="yearly"){
   rebal.names = seq( (firstyear+estyears),lastyear+1,1) 
}else{

   # Name labels using quarters:
   rebal.names = paste(rep( seq( (firstyear+estyears),lastyear,1) , each=4),c("Q1","Q2","Q3","Q4"),sep="")  
   rebal.names = c( rebal.names , paste( lastyear+1, "Q1" , sep="" ) )  
   rebal.names = rebal.names[firstquarter:(length(rebal.names)-4+lastquarter)]

   # Name labels using months: 
   nominalreturns = TRUE;
   if(nominalreturns){ load(file="monthlyR.RData") }else{ load(file="monthlyR_real.RData") } 
   ep = endpoints(monthlyR,on='quarters')
   # select those for estimation period
   ep.start = ep[1:(length(ep)-estyears*4)]+1
   from = time(monthlyR)[ep.start]
   from = seq( as.Date(paste(firstyear,"-01-01",sep="")), as.Date(paste(lastyear-estyears,"-07-01",sep="")), by="3 month") 
   ep.end = ep[(1+estyears*4):length(ep)]
   to = time(monthlyR)[ep.end]
   rebal.names = as.yearmon(to+1)


}


rownames(weightsS1) = rownames(weightsS2) = rownames(weightsS3) = rownames(weightsS4) = rebal.names;
rownames(weightsS5) = rownames(weightsS6) = rownames(weightsS7) = rebal.names; 

rownames(riskcontS1) = rownames(riskcontS2) = rownames(riskcontS3) = rownames(riskcontS4) = rebal.names; 
rownames(riskcontS5) = rownames(riskcontS6) = rownames(riskcontS7) = rebal.names;  


 colorset = gray( seq(0,(cAssets-1),1)/cAssets ) ; 
#due to rounding, the sum of the risk contributions is sometimes 1 + epsilon: avoid this in plot

riskcontS1 = riskcontS1/rowSums(riskcontS1); riskcontS2 = riskcontS2/rowSums(riskcontS2); 
riskcontS3 = riskcontS3/rowSums(riskcontS3); riskcontS4 = riskcontS4/rowSums(riskcontS4); 
riskcontS5 = riskcontS5/rowSums(riskcontS5); riskcontS6 = riskcontS6/rowSums(riskcontS6); 
riskcontS7 = riskcontS7/rowSums(riskcontS7); 

w.names = c( "US bond" , "S&P 500", "EAFE"  , "GSCI" )
 l = 2
mar1 =c(2,l,2,1.1)
mar2 =c(0,l,2,1)
mar3 = c(3,l+1,3,0.1)
mar4 = c(2,l+1,2,0.1)

# Stacked weights plot: 
postscript('stackedweightsriskcont_benchmark.eps') 
   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)

   par(mar=mar3 , cex.main=1)
   chart.StackedBar2(weightsS2,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 )
   chart.StackedBar2(weightsS5,col=colorset,space=0,  main = namelabels[5], 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 )
   chart.StackedBar2(weightsS7,col=colorset,space=0,  main = namelabels[7], 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 )

   par(mar=mar1 , cex.main=1)
   plot.new()
   legend("center",legend=w.names,col=colorset,lwd=6,ncol=4)



   par(mar=mar3 , cex.main=1)
   chart.StackedBar2(riskcontS2,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 )
   chart.StackedBar2(riskcontS5,col=colorset,space=0,  main = namelabels[5], 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 )
   chart.StackedBar2(riskcontS1,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 )

dev.off()

postscript('MinCVaR_alternatives.eps') 
   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)

   par(mar=mar3 , cex.main=1)
   chart.StackedBar2(weightsS3,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 )
   chart.StackedBar2(weightsS4,col=colorset,space=0,  main = namelabels[4], 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 )
   chart.StackedBar2(weightsS6,col=colorset,space=0,  main = namelabels[6], 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 )

   par(mar=mar1 , cex.main=1)
   plot.new()
   legend("center",legend=w.names,col=colorset,lwd=6,ncol=4)
   par(mar=mar3 , cex.main=1)

   chart.StackedBar2(riskcontS3,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 )
   chart.StackedBar2(riskcontS4,col=colorset,space=0,  main = namelabels[4], 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 )
   chart.StackedBar2(riskcontS6,col=colorset,space=0,  main = namelabels[6], 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 )

dev.off()


#############################################################################################
# Plot the CVaR of the portfolios for each period, relatively to CVaR minimum CVaR portfolio


source("R_Allocation/Risk_budget_functions.R"); 
library(zoo); library(fGarch); library("PerformanceAnalytics"); 

# downside risk
alpha = alphariskbudget = 0.05;

#   Load the data

nominalreturns = T;
if(nominalreturns){ load(file="monthlyR.RData") }else{ load(file="monthlyR_real.RData") } 
R = monthlyR

# Define rebalancing periods:

ep = endpoints(monthlyR,on='quarters')
# select those for estimation period
ep.start = ep[1:(length(ep)-estyears*4)]+1
from = time(monthlyR)[ep.start]; from = seq( as.Date(paste(firstyear,"-01-01",sep="")), as.Date(paste(lastyear-estyears,"-07-01",sep="")), by="3 month") 
ep.end = ep[(1+estyears*4):length(ep)]
to = time(monthlyR)[ep.end]
cPeriods = length(from);

# Loop where for each rebalancing period:
# - Compute total CVaR of portfolio

mCVaR = mMU = c()

for( per in c(1:cPeriods) ){
 
 # At the end of each month, we compute the CVaR
 enddates = na.omit(time(window( monthlyR , start = as.Date(to[per]) , end = as.Date(to[per]+90) ),on='months')[1:3])

 for( enddate in enddates ){    

       # add a loop over the next months except when per = cPeriods or -1,-2

       # Estimate GARCH model with data from inception

       inception.R = window(R, start = as.Date(from[1]) , end = enddate );

       # Estimate comoments of innovations with rolling estimation windows
       in.sample.R = window(R, start = as.Date(from[per]) , end = as.Date(to[per]) );
       in.sample.R = checkData(in.sample.R, method="matrix"); 

       # Estimation of mean return
       M = c();
       library(TTR)
       Tmean = 47 # monthly returns: 4 year exponentially weighted moving average
       for( i in 1:cAssets ){
         M = cbind( M , as.vector( EMA(x=inception.R[,i],n=Tmean) ) ) #2/(n+1)
       }
       M = zoo( M , order.by=time(inception.R) )

       # Center returns (shift by one observations since M[t,] is rolling mean t-Tmean+1,...,t; otherwise lookahead bias)
       inception.R.cent = inception.R;
       ZZ = matrix( rep(as.vector( apply( inception.R[1:Tmean, ] , 2 , 'mean' )),Tmean),byrow=T,nrow=Tmean);
       inception.R.cent[1:Tmean,] = inception.R[1:Tmean, ] - ZZ;
       if( nrow(inception.R)>(Tmean+1) ){ 
                 A = M[Tmean:(nrow(inception.R)-1),];
                 A = zoo( A , order.by = time(inception.R[(Tmean+1):nrow(inception.R), ])) ; #shift dates; otherwise zoo poses problem
                 inception.R.cent[(Tmean+1):nrow(inception.R), ] = inception.R[(Tmean+1):nrow(inception.R), ] - A}
       
       # Garch estimation 
       S = c();
       for( i in 1:cAssets ){
            gout =  garchFit(formula ~ garch(1,1), data = inception.R.cent[,i],include.mean = F, cond.dist="QMLE", trace = FALSE )
            if( as.vector(gout@fit$coef["alpha1"]) < 0.01 ){
               sigmat = rep( sd( as.vector(inception.R.cent[,i])), length(inception.R.cent[,i]) ); 
            }else{
               sigmat = gout@sigma.t
            }
            S = cbind( S , sigmat)
       }
       S = zoo( S , order.by=time(inception.R.cent) )

       # Estimate correlation, coskewness and cokurtosis matrix locally using cleaned innovation series in three year estimation window
       selectU = window(inception.R.cent, start = as.Date(from[per]) , end = as.Date(to[per]) )
       selectU = selectU/window(S, start = as.Date(from[per]) , end = as.Date(to[per]) );
       selectU = clean.boudt2(selectU , alpha = 0.05 )[[1]];
       Rcor = cor(selectU)
       D = diag( as.vector(tail(S,n=1)  ),ncol=cAssets )
       sigma = D%*%Rcor%*%D

       # we only need mean and conditional covariance matrix of last observation
       mu = matrix(tail(M,n=1),ncol=1 ) ;
       D = diag( as.vector(as.vector(tail(S,n=1) ) ),ncol=cAssets )
       sigma = D%*%Rcor%*%D
       in.sample.T = nrow(selectU);
       # set volatility of all U to last observation, such that cov(rescaled U)=sigma 
       selectU = selectU*matrix( rep(as.vector(tail(S,n=1)),in.sample.T  ) , ncol = cAssets , byrow = T )
       M3 = PerformanceAnalytics:::M3.MM(selectU)
       M4 = PerformanceAnalytics:::M4.MM(selectU)

       CVaR_period = c( operPortMES(as.numeric(weightsS1[per,]),mu=mu,alpha=alphariskbudget,sigma=sigma,M3=M3,M4=M4)[[1]] ,
          operPortMES(as.numeric(weightsS2[per,]),mu=mu,alpha=alphariskbudget,sigma=sigma,M3=M3,M4=M4)[[1]] ,
          operPortMES(as.numeric(weightsS3[per,]),mu=mu,alpha=alphariskbudget,sigma=sigma,M3=M3,M4=M4)[[1]] ,
          operPortMES(as.numeric(weightsS4[per,]),mu=mu,alpha=alphariskbudget,sigma=sigma,M3=M3,M4=M4)[[1]] ,
          operPortMES(as.numeric(weightsS5[per,]),mu=mu,alpha=alphariskbudget,sigma=sigma,M3=M3,M4=M4)[[1]] ,
          operPortMES(as.numeric(weightsS6[per,]),mu=mu,alpha=alphariskbudget,sigma=sigma,M3=M3,M4=M4)[[1]] ,
          operPortMES(as.numeric(weightsS7[per,]),mu=mu,alpha=alphariskbudget,sigma=sigma,M3=M3,M4=M4)[[1]]  )

       mu_period = c( sum(as.numeric(weightsS1[per,])*mu) , sum(as.numeric(weightsS2[per,])*mu) , 
                      sum(as.numeric(weightsS3[per,])*mu) , sum(as.numeric(weightsS4[per,])*mu) ,
                      sum(as.numeric(weightsS5[per,])*mu) , sum(as.numeric(weightsS6[per,])*mu) ,
                      sum(as.numeric(weightsS7[per,])*mu)  )

       mCVaR = rbind( mCVaR , CVaR_period )
       mMU   = rbind( mMU   , mu_period )
   }
}

colnames(mCVaR) = colnames(mMU) = names[1:7]
mCVaR = xts( mCVaR , 
  order.by = as.Date(time( window( monthlyR , start = as.Date(to[1]) , end = as.Date(to[cPeriods]) ))+1))
mMU = xts( mMU , 
  order.by = as.Date(time( window( monthlyR , start = as.Date(to[1]) , end = as.Date(to[cPeriods]) ))+1))

head(mCVaR[,1:7],2)
#> head(mCVaR[,1:7],2)
#           EqualWeight MinRisk MinRisk_PositionLimit MinRisk_RiskLimit MinRiskConc MinRiskConc_PositionLimit EqualRisk
#1984-01-01      0.0492  0.0320                0.0380            0.0333      0.0352                    0.0388    0.0352
#1984-02-01      0.0490  0.0304                0.0374            0.0322      0.0344                    0.0383    0.0344
save(mCVaR, file="mCVaR.Rdata")


head(mMU[,1:7],2)
#> head(mMU[,1:7],2)
#           EqualWeight    MinRisk MinRisk_PositionLimit MinRisk_RiskLimit MinRiskConc MinRiskConc_PositionLimit   EqualRisk
#1984-01-01 0.009865862 0.00959764            0.01084905       0.009814883 0.009944563                0.01012198 0.009945766
#1984-02-01 0.010176409 0.01003861            0.01105078       0.010193260 0.010297694                0.01042968 0.010299787
save(mMU, file="mMU.Rdata")

###################################################

load(file="mCVaR.Rdata")
load(file="mMU.Rdata")

postscript(file="portfolioMeanCVaR.eps")
par(mfrow=c(2,1),mar=c(3,2,3,2))



plot( mMU[,1]*12  ,  type = "l" , ylim=c(min(mMU),max(mMU))*12,col="darkgray", lwd=1.5 , 
   main = "Expected annualized portfolio return" )
lines( mMU[,2]*12  , type = "l", col="black",lwd=2 , lty=3)
lines( mMU[,7]*12  , type = "l", col="darkgray",lwd=4)
lines( mMU[,5]*12  , type = "l", col="black", lwd=1.5)

legend("bottomleft", legend = c("Equal-Weight","Min CVaR Concentration","Min CVaR+ERC constraint","Min CVaR" ), 
   col=c("darkgray","black","darkgray","black"), lty=c("solid","solid","solid","dashed"), lwd=c(2,2,4,2)  ,cex=0.7) 

plot( mCVaR[,1]  ,  type = "l" , ylim=c(0,max(mCVaR)),col="darkgray", lwd=1.5 , main = "Portfolio CVaR" )
lines( mCVaR[,2]  , type = "l", col="black",lwd=1.5 , lty=3)
lines( mCVaR[,7]  , type = "l", col="darkgray",lwd=4)
lines( mCVaR[,5]  , type = "l", col="black", lwd=1.5)
dev.off()

# do not plot the last month such that it is fully comparable with out-of-sample plots

sel = c( 1 : (nrow(mCVaR)-1) );


# Bear periods
sp500 = window (monthlyR , start = from[1] , end = to[ length(to) ] )[,2]
bear = c(1:length(sp500))[sp500<mean(sp500)]
bear = c(1:length(sp500))[sp500<(-0.12)]
m.bear.dates = list();
i=1;
for( b in bear){
 m.bear.dates[[i]] = c( b-0.5, b+0.5)
 i = i + 1;  
}

out = table.Drawdowns(sp500,top=10) 
start.bear = out$From[out$Depth<(-0.12)]
end.bear = out$Trough[out$Depth<(-0.12)]
start.bear.index = c(1:length(sp500))[ time(sp500) ]
m.bear.dates = list()
v.bear.dates = c()
for( i in 1:length(start.bear) ){
   m.bear.dates[[i]] = c( as.yearmon(start.bear[i]) , as.yearmon(end.bear[i]) )
   v.bear.dates = c( v.bear.dates , seq(start.bear[i],end.bear[i],"days") )
}
v.bear.dates = as.Date( v.bear.dates )


postscript(file="portfolioCVaR.eps")
par(mfrow=c(2,1),mar=c(3,4,1,2))
 
chart.TimeSeries( mCVaR[sel,c(1,5,2)]  , ylim=c(0,max(mCVaR)), ylab = "Portfolio CVaR" , main = "",
 col=c("black","black","darkgray"), lty=c("dashed","solid","solid"), lwd=c(2,2,2,2) , 
   auto.grid = TRUE, minor.ticks = FALSE ,
   period.areas = m.bear.dates , period.color="lightgray",
  date.format.in = "%Y-%m-%d",date.format = "%b %Y") 

legend("topleft", legend = namelabels[c(1,5,2)], 
   col=c("black","black","darkgray"), lty=c("dashed","solid","solid"), lwd=c(2,2,2)  ,cex=0.7) 

chart.TimeSeries( mCVaR[sel,c(4,3,6)]  ,  type = "l" , ylim=c(0,max(mCVaR)), ylab = "Portfolio CVaR"  , main = "",
 col=c("black","black","darkgray"), lty=c("dashed","solid","solid"), lwd=c(2,2,2,2)  ,
   auto.grid = TRUE, minor.ticks = FALSE ,
   period.areas = m.bear.dates , period.color="lightgray",
  date.format.in = "%Y-%m-%d",date.format = "%b %Y") 

legend("topleft", legend = namelabels[c(4,3,6)], 
   col=c("black","black","darkgray"), lty=c("dashed","solid","solid"), lwd=c(2,2,2)  ,cex=0.7) 


dev.off()