Value-at-Risk Historical Bootstrap

If any issues, questions or suggestions feel free to reach me out via e-mail wieczynskipawel@gmail.com or Linkedin. You can also visit my Github.

if(!require(pacman)) install.packages('pacman')
pacman::p_load(PolishStock, dplyr, tidyr, ggplot2, doParallel)

In that blog post we explained Value-at-Risk concept and simulated 10000 portfolios of 4 iShares ETF’s (EEM, IAU, URTH, TIPS) in order to select the optimal ones. However, we calculated VaR and cVaR only one day ahead. In this post we will calculate these measures \(n\) days ahead using bootstrap samples from historical returns.

We use PolishStock package to download data from the website stooq.pl. Daily data for period January 2012 - February 2022 consists of 2345 observations.

tickers <- c('EEM.US', 'IAU.US', 'URTH.US', 'TIP.US')

# for (ticker in tickers){
#   stooq_download(ticker, '20120101', '20220228'
#                  ,destination = paste0(getwd(), '/datasets/'))
# }

close_prices <- df_prices(paste0(tickers, '.csv')
                          ,source = 'datasets/')

returns <- sapply(select(close_prices, -Date)
                  ,FUN = function(x) diff(log(x), lag = 1))

We will simulate \(10000\) portfolios with random weights from Uniform[0,1] distribution. Confidence level for Value-at-Risk is \(\alpha = 95\%\), risk-free interest rate (for generalized Sharpe ratio) is \(\mu_F = 2\%\). For every random portfolio we take \(10000\) samples \(21\) days ahead (approx. 1 month ahead in terms of trading days). To speed up our calculations we use doParallel package in order to utilize multiple CPU cores.

num_of_portfolios <- 1e4
bootstrap_size <- 1e4
days_ahead <- 21
alpha <- 0.95
rf <- 0.02

my.cluster <- makeCluster(detectCores() - 1, type = 'PSOCK')
registerDoParallel(cl = my.cluster)

set.seed(2137)
ahead_metrics <- foreach(
  i = 1:num_of_portfolios
  ,.packages = 'dplyr'
  ) %dopar% {

  portfolio_metrics <- data.frame(
    day = 1:days_ahead
    ,Return = rep(NaN, days_ahead)
    ,VaR = rep(NaN, days_ahead)
    ,VaR_Sharpe = rep(NaN, days_ahead)
    ,cVaR = rep(NaN, days_ahead)
    ,cVaR_Sharpe = rep(NaN, days_ahead)
  )
  
  weights <- c()
  random_weights <- runif(4)
  random_weights <- random_weights / sum(random_weights)
  weights <- random_weights
  names(weights) <- colnames(returns)
  
  portfolio_returns <- returns %*% random_weights
  
  bootstrap_returns <- sample(portfolio_returns
                              ,size = days_ahead * bootstrap_size
                              ,replace = TRUE) %>%
  matrix(nrow = days_ahead, ncol = bootstrap_size) %>%
  apply(2, cumsum)
  
  for (day in 1:days_ahead) {
    portfolio_metrics$Return[day] <- (mean(bootstrap_returns[day,]) + 1)^(252/day) - 1
    portfolio_metrics$VaR[day] <- -quantile(bootstrap_returns[day,], 1 - alpha)
    portfolio_metrics$VaR_Sharpe[day] <- (portfolio_metrics$Return[day] - rf) / portfolio_metrics$VaR[day]
    cVaR_indicator <- (portfolio_returns < -portfolio_metrics$VaR[day])
    portfolio_metrics$cVaR[day] <- -sum(portfolio_returns * cVaR_indicator) / sum(cVaR_indicator)
    portfolio_metrics$cVaR_Sharpe[day] <- (portfolio_metrics$Return[day] - rf) / portfolio_metrics$cVaR[day]
  }
  
  list(weights, portfolio_metrics)
  
}

stopCluster(cl = my.cluster)

Unlisting the results for data vizualization purposes.

d <- days_ahead
d_metrics <- data.frame()
for (i in 1:num_of_portfolios) {
  d_metrics <- bind_rows(
    d_metrics, c(ahead_metrics[[i]][[1]], ahead_metrics[[i]][[2]][d, -1]) %>%
      unlist()
  )
}

quantile_metrics <- data.frame(
  day = 1:days_ahead
  ,Return_q_25 = rep(NA, days_ahead)
  ,Return_q_50 = rep(NA, days_ahead)
  ,Return_q_75 = rep(NA, days_ahead)
  ,VaR_q_25 = rep(NA, days_ahead)
  ,VaR_q_50 = rep(NA, days_ahead)
  ,VaR_q_75 = rep(NA, days_ahead)
  ,cVaR_q_25 = rep(NA, days_ahead)
  ,cVaR_q_50 = rep(NA, days_ahead)
  ,cVaR_q_75 = rep(NA, days_ahead)
)

for (d in 1:days_ahead) {
  temp <- data.frame()
  for (i in 1:num_of_portfolios) {
    temp <- bind_rows(temp, ahead_metrics[[i]][[2]][d, c(2, 3, 5)])
  }
  quantile_metrics[d, 2:4] <- quantile(temp$Return, probs = c(0.25, 0.50, 0.75), na.rm = TRUE)
  quantile_metrics[d, 5:7] <- quantile(temp$VaR, probs = c(0.25, 0.50, 0.75), na.rm = TRUE)
  quantile_metrics[d, 8:10] <- quantile(temp$cVaR, probs = c(0.25, 0.50, 0.75), na.rm = TRUE)
}

In the below chart we see how Expected Return (green) VaR (orange) and cVaR (red) changes with time horizon. Ribbons are interquantile ranges of \(10000\) random portfolios after bootstrap, while the centers of the ribbons are medians of the respective measures. It’s obvious than risk increases with time horizon. Morever, uncertainty connected with risk (ribbon width) also increases with time.

ggplot(quantile_metrics) +
  geom_line(aes(x = day, y = Return_q_50), color = 'green', size = 1) +
  geom_ribbon(aes(x = day, ymin = Return_q_25, ymax = Return_q_75)
              ,fill = 'green' ,alpha= 0.10) +
  geom_line(aes(x = day, y = VaR_q_50), color = 'orange', size = 1) +
  geom_ribbon(aes(x = day, ymin = VaR_q_25, ymax = VaR_q_75)
              ,fill = 'orange' ,alpha= 0.25) +
  geom_line(aes(x = day, y = cVaR_q_50), color = 'red', size = 1) +
  geom_ribbon(aes(x = day, ymin = cVaR_q_25, ymax = cVaR_q_75)
              ,fill = 'red' ,alpha= 0.25) +
  theme_bw() +
  labs(y ='')

Below we see optimal portfolios:
- Minimal VaR and cVaR portfolios have only about \(2.7\%\) - \(3.5\%\) expected annual return and we are 95% confident that monthly (21-sessions) loss won’t exceed \(2.2\%\). In worst case scenario expected monthly dropdown is about \(2.33\%\) - \(2.73\%\). These portfolios mainly consist of treasury inflation-protected bonds.
- Maximal sharpe portfolios have \(8.4\%\) expected annual return and almost two times higher risk than minimal risk portfolios. These portfolios consist of developed markets stocks (\(54\%\)), treasury inflation-protected bonds (\(30\%\)) and gold (\(16\%\)).

optimal_portfolios <- data.frame(
  Min_VaR = round(100*d_metrics[which.min(d_metrics$VaR), ], 2) %>% as.numeric()
  ,Max_VaR_Sharpe = round(100*d_metrics[which.max(d_metrics$VaR_Sharpe), ], 2) %>% as.numeric()
  ,Min_cVaR = round(100*d_metrics[which.min(d_metrics$cVaR), ], 2) %>% as.numeric()
  ,Max_cVaR_Sharpe = round(100*d_metrics[which.max(d_metrics$cVaR_Sharpe), ], 2) %>% as.numeric()
  ) %>% t() %>% as.data.frame()

colnames(optimal_portfolios) <- colnames(d_metrics)
optimal_portfolios

Below we see simulated portfolios on 2-dimensional risk-return planes. Minimal risk portfolios are green triangles and maximal sharpe portfolios are the orange ones.

p1 <- ggplot(d_metrics, aes(x = VaR, y = Return, color = VaR_Sharpe)) +
  geom_point() +
  geom_point(aes(x = (optimal_portfolios$VaR[1] / 100)
                 ,y = (optimal_portfolios$Return)[1] / 100)
             ,color = '#2ca02c', size = 5, shape = 17) +
  geom_point(aes(x = (optimal_portfolios$VaR[2] / 100)
                 ,y = (optimal_portfolios$Return)[2] / 100)
             ,color = '#ff7f0e', size = 5, shape = 17) + 
  theme_bw()

p2 <- ggplot(d_metrics, aes(x = cVaR, y = Return, color = cVaR_Sharpe)) +
  geom_point() +
  geom_point(aes(x = (optimal_portfolios$cVaR[3] / 100)
                 ,y = (optimal_portfolios$Return)[3] / 100)
             ,color = '#2ca02c', size = 5, shape = 17) +
  geom_point(aes(x = (optimal_portfolios$cVaR[4] / 100)
                 ,y = (optimal_portfolios$Return)[4] / 100)
             ,color = '#ff7f0e', size = 5, shape = 17) + 
  theme_bw()

gridExtra::grid.arrange(p1, p2)