This documentation describes the procedure of simulating ETF Prices, for the purpose of finding the optimal weights for risk budgeting Portfolios. The simulation is carried out under strict probabilistic assumption that follows DCC-GJR-GARCH(1,1) Model, that is only relevant to risk parity portfolios, a special type of risk budgeting portfolio. However, alternative investment portfolio strategies will be compared and contrasted with the risk parity portfolio, in the latter part of this documentation. Readers are recommended to refer to the original paper.

1. Importing Packages

The following packages are imported.

# Import packages
rm(list = ls())
if (!require("pacman")) install.packages("pacman")
pacman::p_load(rmgarch, quantmod, abind, doSNOW, quadprog, riskParityPortfolio, CVXR, nlshrink, tidyr, ggridges, ggplot2, reticulate, xts, dplyr, RColorBrewer, formattable, xtable, PerformanceAnalytics)

2. Data

We import the data of historical daily adjusted closing prices data of 11 different exchange-traded funds (ETFs) from Yahoo! Finance. The R package BatchGetSymbols was used for the tickers representing each of the ETF (1st column in the table below). The dates on the price data range from 01-04-2000 to 12-30-2020 (total of 4333 trading days). The asset prices and the dates are concatenated into a xts object suitable for appropriate time series modeling.

# Read in the data
asset_info = read.csv("C:/Users/swcho/Documents/Complete_Example/etfs_countries/asset_info.csv")
dates = read.csv("C:/Users/swcho/Documents/Complete_Example/etfs_countries/etf_dates.csv", header = F)
dates = as.Date(unlist(dates), origin = "1970-01-01") - 719529
prices = read.csv("C:/Users/swcho/Documents/Complete_Example/etfs_countries/etf_countries_prices.csv", header = F)

# Turn into an xts object suitable for future time series modeling
prices = xts(prices, dates)
colnames(prices) = asset_info$Name
n_assets = ncol(prices)

The general overview the asset information containing tickers, the name of the ETFs, and their respective country of origin is summarized as a chart below.

asset_info
##       Ticker                         Name Country
## 1       ^NDX                   NASDAQ 100     USA
## 2      ^GSPC                      S&P 500     USA
## 3       ^HSI              Hang Seng Index     HKG
## 4      ^FTSE                     FTSE 100     GBR
## 5       ^DJI Dow Jones Industrial Average     USA
## 6     ^GDAXI        DAX Performance-Index     GER
## 7       ^RUT                 Russell 2000     GBR
## 8      ^FCHI                       CAC 40     FRA
## 9      ^BVSP                     Ibovespa     BRA
## 10 000001.SS          SSE Composite Index     CHN
## 11     ^N225                   Nikkei 225     JPN

We then plot a line graph of the evolution of historical adjusted closing prices for 11 assets. The dates and the historical asset prices go on the x and y axis, respectively.

# Data Visualization

p1 = plot(xts(prices, dates)
            , legend.loc = "topleft"
            , main="Historical Asset Prices")

print(p1)

3. Market Model

Before diving into the details of the model framework, We denote log returns as \(R_{i,t}=log\left(\frac{P_{i,t}}{p_{i,t-1}}\right)\) for \(i = 1,...,11\) and \(t = 2,...,4333\). The historical asset prices are manipulated to create a new object, ‘returns’ containing historical log returns of each ETF. Note that the very first row is NA.

# NA in the first row
returns = diff(log(prices))
data.frame(head(returns, 5))

For each of the 11 ETFs, a univariate GARCH model is fitted to model the heteroskedasticity of the (log) returns. Some Probabilistic assumptions made for each of the univariate GARCH model are as follows:

The marginal univariate GARCH models are then merged together fit a DCC-GJR-GARCH(1,1) model, which models the joint distribution of the returns among the 11 assets. DCC-GJR-GARCH is chosen as it gives stronger weights to negative shocks than positive ones at time t-1. The conditional variance of the returns at time t is modeled as follows:

\[h_{ii,t} = \theta_i+(a_i+c_i[\epsilon_{i,t-1}<0])\epsilon^2_{i,t-1} + b_ih_{ii,t-1}\] * The square brackets, [] represents an indicator function. * \(a_i, b_i, c_i\geq0\) and \(\epsilon_{i,t}\) are i.i.d standard normally distributed.

3.1. Parameter Specifications

We select the window of model fitting as 1260 days or 5 years (assuming 252 trading days per year). The goal is to simulate the price 5 days ahead. The simulation is run 1000 times for each of the window (1000 simulations for a simulated day). The reason for simulating the prices is that we use those simulated prices, say at time t+5, to calculate the expected shortfall (over the horizon of 5 days) at a selected singfinicance level (\(\alpha = 0.95\), unless explicitly stated otherwise). For the sake of simplicity, we obtain predictions for the first 100 days.

n_fit = 5*252         # number of data points used to fit the models (1260)
n_days_ahead = 5      # number of days ahead for computing the risk measure 
n_simulations = 1000  # number of simulations

n_start = 100         # "burn-in" for the GARCH estimation


#subset the data
prices <- prices[c(1:(n_fit+100)), ]
dates <- dates[c(1:(n_fit+100))]

# NA in the first row
returns = diff(log(prices))
data.frame(head(returns, 5))
##              NASDAQ.100       S.P.500 Hang.Seng.Index     FTSE.100
## 2000-01-04           NA            NA              NA           NA
## 2000-01-05 -0.011027215  0.0019203377     -0.07452522 -0.019694921
## 2000-01-06 -0.048636068  0.0009552219     -0.04474878 -0.013664084
## 2000-01-07  0.054971268  0.0267299501      0.01651928  0.008894378
## 2000-01-11  0.004170236 -0.0020207518      0.02919957  0.002165299
##            Dow.Jones.Industrial.Average DAX.Performance.Index Russell.2000
## 2000-01-04                           NA                    NA           NA
## 2000-01-05                 0.0112765600           -0.01296989  0.000940195
## 2000-01-06                 0.0116742387           -0.00418432 -0.007315272
## 2000-01-07                 0.0236489718            0.04618244  0.026920117
## 2000-01-11                -0.0009967589            0.01613381  0.008767310
##                  CAC.40     Ibovespa SSE.Composite.Index   Nikkei.225
## 2000-01-04           NA           NA                  NA           NA
## 2000-01-05 -0.034494928  0.024552579         0.002351544 -0.024521328
## 2000-01-06 -0.005414624 -0.008531209         0.037768657 -0.020391493
## 2000-01-07  0.016288310  0.012463143         0.035340819  0.001382809
## 2000-01-11  0.013321571  0.016057762        -0.024579519  0.035502266
n_days = nrow(prices)  # Only performed 100 times for the sake of simplicity
step_seq = seq(n_fit+1, n_fit + 100)
n_steps = length(step_seq)

Next, parameter specifications for univariate GARCH(1,1) and DCC-GJR GARCH(1,1) are made.

# Mean Model Specifications (garch orders)

ar = 0 #Auto Regressive
ma = 0 # Moving Average

#Variance Model Specifications
garch_model = "gjrGARCH"
p_garch = 1
q_garch = 1

# Conditional density to use for the innovations
marginal_dist = "norm" #"std"

# DCC Specifications
multi_dist = "mvnorm" #"mvt"
p_dcc = 1
q_dcc = 1

# Garch solver
garch_solver_1 = "gosolnp"

# The “gosolnp” solver allows for the initialization of multiple restarts of the solnp solver with randomly generated parameters

# Model Specification

# Marginal GARCH models (p,q) for all assets
unispec = ugarchspec(mean.model = list(armaOrder = c(ar,ma)),
                     variance.model = list(model = garch_model,
                                           garchOrder = c(p_garch, q_garch)),
                     distribution.model = marginal_dist) 

garch_spec = multispec(replicate(n_assets, unispec)) 

# garch_spec replicates this specification for all assets
# Each component will have its own parameters

# DCC model
spec = dccspec(uspec = garch_spec, 
               dccOrder = c(p_dcc, q_dcc), 
               distribution = multi_dist)

# DCC Simulation random seed
# note that we need to use the random seed in the dccsim function, not in set.seed.
seed = 1

3.2. Model Fitting and Simulations

We create an empty list, in order to stack the 100 days worth of data of 1000 price simulations for 11 assets.

# First 100 Days =====================
stack <-  list() # Stacking Data Arrays

garch_solver = garch_solver_1

# Finding the parameters for marginal GARCH,  using the chosen algorithm via Maximum Likelihood Estimation

We then create a function for fitting a model for a single asset, for a selected window.

compute_fit = function(j){
  returns_window_j = returns_window[, j] # reminder: the first row of the returns is NA
  prices_window_j = prices_window[, j]
  
  unifit = ugarchfit(unispec, returns_window_j
                       , solver = garch_solver
                       , solver.control = list(rseed = seed
                                             # , pars=fitlist[[j]]@fit$coef * 1000
                                             # , pars = c(0.0003, 0.000008608, 0.04, 0.9)/10
                                             # , trace = T
                       ))
  
  return(unifit)}

The algorithm of the for-loop is as follows:

  1. Select the first window of 1260 days (starting with the very first window that spans from Day 2 to Day 1261) of the the returns (object with historical log returns) and prices (object with historical ETF prices).

  2. Estimate the parameters and fit univariate GARCH(1,1) for 11 assets within the selected window, resulting in 11 model fit objects.

  3. Estimate the parameters and fit a DCC-GJR GARCH(1,1) using the 11 fitted univariate GARCH(1,1) models. This is the model that is used to simulate.

  4. Simulate (1000 times) the log returns for the five days ahead of Day 1261 (Day 1262 to Day 1266). We should end up with 11 data arrays with the dimension of 1000 simulations x 5 days.

  5. Add the 5 days worth of simulated log returns, to get the cumulative log returns over the span of 5 days (Use of the logarithmic rule, \(ln(AB) = ln(A) + ln(B)\)).

  6. Apply exponential function across the rolled up log returns to get the 5 day compounded returns. For each of the 11 assets, we should have 1000 simulations of 5 day compounded returns.

  7. Multiply the 5 day compounded return values to the final historical price value within the selected window to produce a prediction of price 5 days ahead. For each of the 11 assets, we should now have 1000 simulations of 5-day-ahead price predictions (For the first window, these are price predictions for day 1266). Stack the prediction results in a list.

  8. Shift the window for model fitting by a day. Iterate 100 times.

for(step in step_seq){
  index_in_the_window = (step-n_fit+1):(step) 
  # Very first iteration should be rows 2 to 1261 of the returns xts object
  returns_window = returns[index_in_the_window, ] # reminder: the first row of the returns is NA
  prices_window = prices[index_in_the_window, ]

  fitlist = lapply(1:n_assets, compute_fit) # Apply compute_fit across 11 assets in a selected window
  desc = list()
  desc$type = "equal"
  desc$asset.names = colnames(returns_window)
  
  multf = new("uGARCHmultifit",
             fit  = fitlist,
             desc = desc)

  #  Produce Parameters for the joint distribution
  fit = dccfit(spec, 
             data = returns_window, 
             solver = garch_solver,
             fit.control = list(eval.se = FALSE) 
             , solver.control = list(rseed = seed)    
             , fit = multf)


  # Simulation ========================
  dcc_sim_out = dccsim(fitORspec = fit, n.sim = n_days_ahead, n.start = n_start, 
                       m.sim = n_simulations, rseed = seed )
  
  # Simulate the price 5 days ahead, 1000 times using the parameters estimated
  # Using the dcc model built, we simulate returns for 5 days ahead

  returns_sim = abind(dcc_sim_out@msim$simX, along = 3)
  returns_sim = aperm(returns_sim, c(3,1,2)) # 1000 simulations x 5 days x 11 assets
  returns_sim_end = apply(returns_sim, c(1,3), sum)
  # 1000 simulated, 5 day compounded return, for each of the 11 assets
  
  
  prices_sim = sweep( 
    x = exp(returns_sim_end), # exp of the cumulative return in n_days_ahead (dim = n_simulations x n_assets)
    MARGIN = 2,
    FUN = "*",
    STAT = prices_window[n_fit, ] # length(index_in_the_window) = n_fit
  )
    
  prices_sim = array(prices_sim, dim = c(1, dim(prices_sim)))
 
  # Stack the results in the empty list created in the beginning
  
  i = step - n_fit
  
  stack[[i]] <- prices_sim
  
  print(paste0('Progress ',(i/n_steps*100),'%'))
  
}

# End of loop

By the end of the loop, we have a list object, stack, containing 100 days worth of predictions (in order), each day with 1000 simulations for 11 assets. We can visualize the results from the simulations.

# Let's visualize the very first simulation out of total of 1000 simulations

sim_val <- rep(0, (n_assets*n_steps))
sim_val <- matrix(sim_val, ncol = 11)

for (i in seq(1:n_steps)){
  sim_val[i, ] <- matrix(stack[[i]], ncol = 11)[1,]
}

dates = read.csv("C:/Users/swcho/Documents/Complete_Example/etfs_countries/etf_dates.csv", header = F)
dates = as.Date(unlist(dates), origin = "1970-01-01") - 719529
dates_sim <- dates[-c(1:(n_fit + n_days_ahead))] 
dates_sim <- dates_sim[c(1:n_steps)] # We have dates from Day 1266 to Day 1365 (100 Days)

sim_val_xts <- xts(sim_val, dates_sim)
colnames(sim_val_xts) <- asset_info$Name
p2 = plot(sim_val_xts
            , legend.loc = "topleft"
            , main= "Simulated Asset Prices (Sim. #1)")

print(p2)

4. Alternate Portfolio Strategies

Risk budgeting portfolios aren’t just diversified in terms of how each asset is allocated among different assets, but they are diversified in terms how the investor wants each asset to contribute to the aggregate portfolio risk. Each asset within the risk budgeting portfolio contributes a prespecified amount to the aggregate portfolio risk. If the predetermined amount of each asset’s contribution is equal across all assets within the portfolio, the portfolio is called “risk parity portfolio.” The simulated price data above, using DCC-GJR GARCH(1,1) is inputted, to create and calculate the optimal weights of a risk parity portfolio (rpb). Within the paper, the performance of the risk parity portfolio is compared to other investment strategies. The details for the assumptions for each of the portfolio is explained in detail in Section 4 of the paper. Alternate investment strategies mentioned in the paper are the following:

The summary of the competing portfolios and the mehtods of obtaining the optimal weights for each of the portfolio is summarized below.

Investment Strategies Abbreviations Optimizing Functions
Maxmmum Sharpe Ratio msr \(\underset{w}{max}(\frac{w^T\mu}{w^T \Sigma w})\) subject to \(w\geq{0}\) and \(1^T w = 1\)
Markowitz Mean-Variance mmv \(\underset{w}{min}({w^T \Sigma w})\) subject to \(w\geq{0}\) and \(1^T w = 1\)
Equal Weights ew Equal weights are allocated to each of the asset and constant throughout time
Gaussian Risk Parity grp -
Risk parity using DCC-GARCH rpa -
Risk Parity using DCC-GJR-GARCH rpb -

Global minimum variance (mmv) portfolio and the standard deviation risk parity (sdrP) have been omitted here. However, the results and the further details for these portfolios are included in the original paper. The Gaussian risk parity portfolio (grp) assumes that the asset returns follow a multivariate Gaussian distribution. The risk parity portfolio using DCC-GARCH(1,1) (DCC model using simple GARCH) (rpa) assumes that the conditional variances of the returns at time are as follows:

\[h_{ii,t} = \theta_{i} + a_{i}\epsilon^2_{i,t-1} + b_ih_{ii,t-1} \] where \(a_i, b_i, c_i\geq0\) and \(\epsilon_{i,t}\) are i.i.d standard normally distributed. Both grp and rpa portfolios measure risk using expected shortfall at level \(\alpha\), as the coherent risk measure. The optimal weights for Gaussian risk parity portfolio (grp) and risk parity portfolio using DCC-GARCH(1,1) (rpa) are found using algorithm in Section 3 of the paper.

First we define a function, which computes the portfolio exposure and wealth, given daily weights and the initial capital. We make an assumption that the investor has an initial endowment of 100 dollars, and invests entirely in the 11 assets. We further assume that all portfolios are rebalanced daily. The exposure and wealth are defined as follows:

Terms Definitions
Exposure The exact dollar amount that is invested in each of the asset
Wealth The portfolio value of the investor

We first define a function which calculates the exposure and the wealth for each of the portfolio given weights and the initial capital.

compute_exposure_and_wealth = function(weights, initial_capital){
  # the input weights should be normalized to sum up to 1.
  n_days = nrow(weights)
  d = ncol(weights)
  
  wealth = rep(NA, n_days)
  exposure = matrix(NA, n_days, d)
  
  wealth[1] = initial_capital
  exposure[1,] = unlist(weights[1,] * initial_capital)
  
  for(t in 2:n_days){
    wealth[t] = sum(exposure[t-1, ] * prices[t,] / prices[t-1,])
    exposure[t, ] = unlist(weights[t,] * wealth[t])
  }

  return(list(exposure = exposure,
              wealth = xts(wealth, dates[1:n_days])))
}

We further define the functions that calculate the weights for maxmimum sharpe ratio, and markowitz mean-variance portfolio. The weights are calculated by finding the solutions to the corresponding optimization functions stated above. The portfolios are adjusted based on a set of 252*5 = 1260 days previous to the current period, time t. Data in the interval \([t-1260, t-1]\) is used to construct the portfolio at time t. The \(\mu\) and \(\Sigma\) are empirical mean and covariance of returns in a selected window of 1260 days.

# The code has been adapted to use as input the returns and not prices
# this is to make things consistent with our return definition for our risk parity

# tangency portfolio (maximum sharpe ratio)
max_sharpe_ratio <- function(returns, Sigma) {
  # relative returns = (P_t - P_{t-1})/P_{t-1}
  N <- ncol(returns)
  mu <- colMeans(returns)
  if (all(mu <= 1e-8))
    return(rep(0, N))
  Dmat <- 2 * Sigma
  Amat <- diag(N)
  Amat <- cbind(mu, Amat)
  bvec <- c(1, rep(0, N))
  dvec <- rep(0, N)
  res <- solve.QP(Dmat = Dmat, dvec = dvec, Amat = Amat, bvec = bvec, meq = 1)
  x <- res$solution
  w = x/sum(x)
  return(w)
}

# define Markowitz mean-variance portfolio
Markowitz_portfolio_fun <- function(returns, Sigma) {
  # relative returns = (P_t - P_{t-1})/P_{t-1}
  X <- returns  # compute log returns
  mu <- colMeans(X)  # compute mean vector
  # design mean-variance portfolio
  w <- Variable(nrow(Sigma))
  prob <- Problem(Maximize(t(mu) %*% w - 0.5*quad_form(w, Sigma)),
                  constraints = list(w >= 0, sum(w) == 1))
  result <- solve(prob)
  return(as.vector(result$getValue(w)))
}

We next, calculate the weights for the ew, msr, mmv portfolios for the first 100 days.

colnames(prices) = asset_info$Ticker

# Computing the competing portfolios ==============
returns = apply(prices, 2, function(x) x[-1]/x[-length(x)] - 1)
returns = rbind(NA, returns) # to be consisted with the RP portfolio code

n_window = 252*5 # should be consistent with our risk parity approach

  msr_weights =    # maximum sharpe ratio (tangency)
  mmv_weights =    # Markowitz mean-variance portfolio
  matrix(NA, nrow(returns), ncol(returns))


for(t in (n_window+1):nrow(returns)){
  ret = returns[(t-n_window+1):t, ]
  Sigma <- cov(ret)

  msr_weights[t,] = max_sharpe_ratio(ret, Sigma)
  mmv_weights[t,] = Markowitz_portfolio_fun(ret, Sigma)
  cat(t, "\n")
}

msr_weights <- msr_weights[-c(1:n_fit),]
mmv_weights <- mmv_weights[-c(1:n_fit),]

# these vectors should be aligned with dates
portfolios = c("ew", "msr", "mmv")

prices <- as.data.frame(prices[-c(1:n_fit)])
dates <- dates[-c(1:n_fit)]
  
initial_capital <- 100

# the last competing portfolio --------------------
ew_weights = matrix(1/n_assets, n_steps, n_assets)

Using the msr, mmv weights, we calculate the daily portfolio wealth and exposures.

#ew, msr, mmv portfolio wealths/exposure

for(portfolio_names in portfolios){
  assign(paste0(portfolio_names, "_out"),
         compute_exposure_and_wealth(get(paste0(portfolio_names,"_weights")), 
                                     initial_capital)
  )
  
  assign(paste0(portfolio_names, "_wealth"),
         get(paste0(portfolio_names, "_out"))$wealth
  )

  # Exposure per asset ============
  assign(paste0(portfolio_names, "_exposure"),
         get(paste0(portfolio_names, "_out"))$exposure
  )
}

We now move on to calculate the daily portfolio weights for grp, rpa, rpb, ew portfolios. For grp, rpa, and rpb portfolios, we utilize Expected Shortfall (ES) as the coherent risk measure. ES at security level \(\alpha \in (0, 1)\) of a random variable \(X\in\mathbb{L^2}\) is defined by the following:

\[ ES_\alpha(X) := \frac{1}{1-\alpha}\int_\alpha^1 F_X^{-1}(u)du\]

The typical values for the security level used in financial applications are \(\alpha=0.95\) and \(\alpha = 0.99\). We include expected shortfall significance levels with \(\alpha \in \left\{0.80, 0.85 0.90, 0.95, 0.96, 0.97, 0.98, 0.99 \right\}\) in order to further explore the portfolio sensitivity to the significance levels.

garch_model_list = c("sGARCH", "gjrGARCH")

for(garch_model in garch_model_list){
  cat(garch_model, "\n")
  asset_info = read.csv("C:/Users/swcho/Documents/Complete_Example/etfs_countries/asset_info.csv")
  dates = read.csv("C:/Users/swcho/Documents/Complete_Example/etfs_countries/etf_dates.csv", header = F)
  dates = as.Date(unlist(dates), origin = "1970-01-01") - 719529
  prices = read.csv("C:/Users/swcho/Documents/Complete_Example/etfs_countries/etf_prices.csv", header = F)

# Turn into an xts object suitable for future time series modeling
  prices = xts(prices, dates)
  colnames(prices) = asset_info$Name
  n_assets = ncol(prices)
  model = paste0("arma_",ar,ma,"_",
               multi_dist,"_dcc_",p_dcc, q_dcc,"_",
               marginal_dist,"_", garch_model , "_",p_garch,q_garch) #"t_dcc_t_garch" #"gauss_dcc_gauss_garch"
  
  # it is easier if we transform the prices object from xts to data.frame
  prices_xts = prices
  prices = as.data.frame(prices)
  
  source("C:/Users/swcho/Documents/Complete_Example/weights_analysis_countries/set_parameters_plot.R")
  source("C:/Users/swcho/Documents/Complete_Example/weights_analysis_countries/3_create_portfolio_1_1.R")
}

all_portfolios <- all_portfolios[! all_portfolios %in% c('gmv','sdrp', 'msr', 'mmv', 'ew')]

for (port in all_portfolios){
  assign(paste0(port, '_weights'), get(paste0(port, '_weights'))[c(1:100), ]) # truncate the weights for just first 100 days
}

We proceed to calculate wealth and exposures for grp, rpa, rpb portfolios.

asset_info = read.csv("C:/Users/swcho/Documents/Complete_Example/etfs_countries/asset_info.csv")
dates = read.csv("C:/Users/swcho/Documents/Complete_Example/etfs_countries/etf_dates.csv", header = F)
dates = as.Date(unlist(dates), origin = "1970-01-01") - 719529
prices = read.csv("C:/Users/swcho/Documents/Complete_Example/etfs_countries/etf_prices.csv", header = F)
prices = xts(prices, dates)
colnames(prices) = asset_info$Name
n_assets = ncol(prices)

prices <- prices[c(1:(n_fit+100)), ]
prices <- as.data.frame(prices[-c(1:n_fit), ]) # Price data from day 1261 to day 1360
dates <- dates[c(1:(n_fit+100))]
dates <- dates[-c(1:n_fit)]

for(portfolio_names in all_portfolios){
  assign(paste0(portfolio_names, "_out"),
         compute_exposure_and_wealth(get(paste0(portfolio_names,"_weights")), 
                                     initial_capital)
  )
  
  assign(paste0(portfolio_names, "_wealth"),
         get(paste0(portfolio_names, "_out"))$wealth
  )

  # Exposure per asset ============
  assign(paste0(portfolio_names, "_exposure"),
         get(paste0(portfolio_names, "_out"))$exposure
  )
}

5. Results

5.1. Histogram of Returns

Using the exposure and the wealth results for each of the competing portfolio, we can further explore their performance via data visualization techniques. First, we utilize partially overlapping ridge line plots in order to visualize the overall portfolio returns distributions of each of the competing portfolio.

all_portfolios_2 <- all_portfolios_2[!all_portfolios_2 %in% c('gmv', 'sdrp')]

all_wealths = NULL
for(portfolio_names in all_portfolios_2){
  all_wealths = cbind(all_wealths,
                      get(paste0(portfolio_names,"_wealth"))
  )
}

all_wealths = xts(all_wealths, dates)
colnames(all_wealths) = all_portfolios_2

all_wealths_plots = all_wealths

all_returns = apply(all_wealths_plots, 2, diff)/all_wealths_plots[-nrow(all_wealths_plots),]
t_all_returns = data.frame(portfolio = colnames(all_returns),
                           t(all_returns))

returns_long = pivot_longer(data = t_all_returns, 
                            cols = -c("portfolio"), 
                            names_to = "date", 
                            values_to = "returns")

i = 0

for(portfolio_names in portfolio_type_plot){
  i = i + 1
  
  returns_long_plot = returns_long %>% 
    filter(portfolio %in% unlist(portfolio_names))
  returns_long_plot$portfolio = factor(returns_long_plot$portfolio, 
                                       levels = rev(unlist(portfolio_names)))
  
  
  p = ggplot(returns_long_plot, 
             aes(x = returns, y = portfolio, group = portfolio)) +
    geom_density_ridges(scale = 5, size = 0.5, rel_min_height = 0.00, aes(fill = portfolio)) +
    scale_fill_cyclical(values = rev(col_type_plot[[i]])) +
    theme_ridges(font_size = 20) +
    scale_x_continuous(limits = c(-0.025, 0.025), expand = c(0, 0)) + #https://stackoverflow.com/questions/43475764/ggplot-non-finite-values-error
    coord_cartesian(clip = "off")
  
  print(p)
  #dev.off()  
}

As we fix the portfolio strategy for grp, rpa or rpb, and vary the significance level \(\alpha\), it can be seen that, there is almost no difference between them. This suggests some robustness with respect to \(\alpha\).

5.2. Plotting Wealth over Time

Next, the values of each portfolio with initial endowment of 100 dollars, are plotted as a time series plot over 100 days from 2006-02-24 to 2006-08-07. For grp portfolio and rpa, rpb portfolios, we calculated the ES at significance level \(\alpha\) = 0.95, using N = 1000 simulations at time t+h. The following portfolios are selected to be plotted:

  • Maximum sharpe ratio portfolio (msr)
  • Markowitz mean-variance portfolio (mmv)
  • Equal weights portfolio (ew)
  • Gaussian risk parity portfolio with expected shortfall (ES) significance level of \(\alpha\) = 0.95 (grp95)
  • Risk parity portfolio using DCC-GARCH on simple GARCH with expected shortfall (ES) significance level of \(\alpha\) = 0.95 (rpa95)
  • Risk parity portfolio using DCC-GARCH on GJR-GARCH with expected shortfall (ES) significance level of \(\alpha\) = 0.95 (rpb95)
all_portfolios_3 <- all_portfolios_2[all_portfolios_2 %in% c('msr', 'mmv', 'ew', 'grp95', 'rpa95', 'rpb95')]

i_plot = 1
all_wealths = NULL
for(portfolio_names in all_portfolios_3){
  all_wealths = cbind(all_wealths,
                      get(paste0(portfolio_names,"_wealth"))
  )
}

all_wealths = xts(all_wealths, dates)
colnames(all_wealths) = all_portfolios_3

all_wealths_plots = all_wealths

plot.xts(all_wealths_plots
            , col = col_portfolios
            , legend.loc = "bottomleft"
            , main="Portfolio Values (Initial endowment of $100)")

The results show that the risk parity portfolios outperform msr, mmv, and ew portfolios in the short term. The long term graph can be found under Section 5 of the original paper.

5.3. Asset Allocation

We proceed to plot barplots of the weights of all 11 assets, across the selected portfolios, msr, mmv, ew, grp95, rpa95, and rpb95, over time. First, we create the legend, showing the corresponding colors for each of the asset on the allocation plot.

for(portfolio_names in all_portfolios_3){

  portfolio_weights = get(paste0(portfolio_names,"_weights"))
  
  data_normalized_barplot = xts(portfolio_weights, dates[1:nrow(portfolio_weights)])
  
  par(mar=c(5, 1, 1, 1), xpd=TRUE)
  my_barplot = barplot(data_normalized_barplot[,ncol(data_normalized_barplot):1], 
                       col = rev(col_assets), 
                       bty = "L",
                       #yaxt = 'n', xaxt = 'n', 
                       border = rev(col_assets)
                       ,main = paste0(portfolio_names)
                       )

  print(my_barplot)
}

In the following weight allocation graph over time, msr and mmv are the most dissimilar portfolios in terms of their daily rebalanced weights across 11 assets. It can be seen that for mmv, it is invested in only one single asset (Ibovespa) for the first 100 days. The results for the msr portfolio is not as extreme but similar, with the portfolio invested in 3 assets at most at one certain point in time. The weights for the equal weights is self-evident. For the risk parity strategies (grp95, rpa95, rpb95), the weights are more balanced over time, with rpa95 and rpb95 being very similar. The weights under rpb95 however, does seem more robust with less fluctuations, suggesting that in order to take full advantage of the risk parity strategy using ES as the coherent risk measure, a more flexible underlying model is required.

References

Freitas Paulo da Costa, Bernardo, Silvana M. Pesenti, and Rodrigo Targino. “Risk Budgeting Portfolios from Simulations.” Silvana M. and Targino, Rodrigo, Risk Budgeting Portfolios from Simulations (February 18, 2022) (2022)