The critical line algorithm (CLA) was developed by Harry Markowitz to solve the corner portfolios of the efficiency frontier in his famous 1952 Portfolio Selection research. An open source version for Python exists and is outlined in the 2013 paper below. CriticalLineAlgo is a lightweight R implementation with minimal imports. The only outside package is MASS for general (Moore-Penrose) matrix inversion.
This package mostly follows the notation of Bailey and Lopez de Prado 2013 and Niedermayer and Niedermayer 2006.
There are two main functions to interface with: run_cla and calc_target_vol. The run_cla function solves for the corner portfolios and output a list with the stored output. The calc_target_vol function takes the list output from run_cla and finds the weights of an efficient portfolio for a specific level of volatility.
In order to run the CLA a mu_vec of estimated returns and cov_mat of estimated covariance are needed. This introduction reads the mu_vec and cov_mat from csv files stored in Google Drive. The estimates are historical mean and cov estimates from 16 ETF time-series to represent different asset classes. The mu_vec and cov_mat are not predictions of the future. Instead they are used here to demonstrate how this package works.
mu_data <- read.csv(sprintf('https://docs.google.com/uc?id=%s&export=download',
'1vUPvKqJSR-3ZlWTWwf1r2kjguEDqmDM6'),
row.names = NULL)
cov_data <- read.csv(sprintf('https://docs.google.com/uc?id=%s&export=download',
'1g53v1D9mXvfcYKMsquHIj2CF7Lrn6-GX'),
row.names = NULL)
asset_class <- mu_data$ETF
mu_vec <- mu_data$mu
cov_mat <- as.matrix(cov_data[, 2:ncol(cov_data)])
mu_data
#> ETF mu
#> 1 IVV 0.14335496
#> 2 IWM 0.12231009
#> 3 VTV 0.13113963
#> 4 IWF 0.16094526
#> 5 EFA 0.05844601
#> 6 VWO 0.02959126
#> 7 TLT 0.05983401
#> 8 IEF 0.03455211
#> 9 EMB 0.05410923
#> 10 LQD 0.05320704
#> 11 HYG 0.05960452
#> 12 VNQ 0.10600503
#> 13 AMJ 0.01488107
#> 14 IGF 0.07625922
#> 15 GNR 0.01742278
#> 16 GLD 0.01517350run_cla is the interface that runs the algorithm and returns the results. The mu_vec and cov_mat are passed in as the first two arguments.
The $wgt_list in the store output contains the weights of the corner portfolios. These weights can easily be combined and used to calculate the corner portfolios’ expected return and volatility.
corner_wgt <- do.call(cbind, store$wgt_list)
exp_ret <- apply(corner_wgt, 2, function(m, w) {t(w) %*% m}, m = mu_vec)
exp_vol <- apply(corner_wgt, 2, function(covar, w) {sqrt(t(w) %*% covar %*% w)},
covar = cov_mat)
plot(sqrt(diag(cov_mat)), mu_vec, col = 'darkgrey',
ylab = 'Estimated Return', xlab = 'Estimated Volatility')
points(exp_vol, exp_ret, col = 'darkgreen', lwd = 2, type = 'b')
title('Efficiency Frontier Example')
mtext('Corner Portfolios in Green, Individual Assets in Grey')
rownames(corner_wgt) <- asset_class
colnames(corner_wgt) <- paste('Vol ', round(exp_vol * 100, 1))
round(corner_wgt * 100, 1)
#> Vol 15 Vol 15 Vol 8.3 Vol 6.8 Vol 6.1 Vol 5 Vol 4.8 Vol 4.8
#> IVV 0 0 0.0 0.0 0.0 0.0 0.0 0.0
#> IWM 0 0 0.0 0.0 0.0 0.0 0.0 0.0
#> VTV 0 0 0.0 0.0 0.0 2.0 4.4 5.0
#> IWF 100 100 59.3 44.3 39.0 28.6 25.2 24.2
#> EFA 0 0 0.0 0.0 0.0 0.0 0.0 0.0
#> VWO 0 0 0.0 0.0 0.0 0.0 0.0 0.0
#> TLT 0 0 40.7 23.9 14.8 0.0 0.0 0.0
#> IEF 0 0 0.0 0.0 14.8 39.5 42.6 43.3
#> EMB 0 0 0.0 0.0 0.0 0.0 0.0 0.6
#> LQD 0 0 0.0 31.8 31.4 29.9 27.9 26.9
#> HYG 0 0 0.0 0.0 0.0 0.0 0.0 0.0
#> VNQ 0 0 0.0 0.0 0.0 0.0 0.0 0.0
#> AMJ 0 0 0.0 0.0 0.0 0.0 0.0 0.0
#> IGF 0 0 0.0 0.0 0.0 0.0 0.0 0.0
#> GNR 0 0 0.0 0.0 0.0 0.0 0.0 0.0
#> GLD 0 0 0.0 0.0 0.0 0.0 0.0 0.0
#> Vol 4.3 Vol 3.9 Vol 3.9 Vol 3.8
#> IVV 0.0 5.5 5.7 0.0
#> IWM 0.0 0.0 0.0 0.0
#> VTV 7.2 7.0 7.0 9.3
#> IWF 14.1 0.9 0.0 0.0
#> EFA 0.0 0.0 0.0 0.0
#> VWO 0.0 0.0 0.0 0.0
#> TLT 0.0 0.0 0.0 0.0
#> IEF 51.1 59.6 59.5 58.7
#> EMB 2.4 4.1 4.1 4.3
#> LQD 13.9 0.0 0.0 0.0
#> HYG 11.3 22.9 23.6 27.6
#> VNQ 0.0 0.0 0.0 0.0
#> AMJ 0.0 0.0 0.0 0.0
#> IGF 0.0 0.0 0.0 0.0
#> GNR 0.0 0.0 0.0 0.0
#> GLD 0.0 0.0 0.0 0.0run_cla supports adding lower and upper bound weight limits. This version of the CLA has a long-only constraint (i.e., the sum of the portfolio weights = 1). There are two ways of entering asset constraints. A single numeric value will be applied to each asset, e.g., entering 0.5 as the upper bound constraint will limit every asset to a maximum weight of 50%. Alternatively, a vector of weights that correspond to each asset in the mu_vec and cov_mat can be entered as input to allow for different constraints for each asset. The code below is an example of the latter: a constraint of 50% will be applied to all assets except HYG, AMJ, and GLD will be constrained at 10%.
up_bound <- rep(0.5, length(mu_vec))
up_bound[asset_class %in% c('HYG', 'AMJ', 'GLD')] <- 0.1
store_constr <- run_cla(mu_vec, cov_mat, up_bound = up_bound)
corner_wgt <- do.call(cbind, store_constr$wgt_list)
exp_ret <- apply(corner_wgt, 2, function(m, w) {t(w) %*% m}, m = mu_vec)
exp_vol <- apply(corner_wgt, 2, function(covar, w) {sqrt(t(w) %*% covar %*% w)},
covar = cov_mat)
plot(sqrt(diag(cov_mat)), mu_vec, col = 'darkgrey',
ylab = 'Estimated Return', xlab = 'Estimated Volatility')
points(exp_vol, exp_ret, col = 'darkgreen', lwd = 2, type = 'b')
title('Efficiency Frontier Example with Constraints')
mtext('Corner Portfolios in Green, Individual Assets in Grey')
rownames(corner_wgt) <- asset_class
colnames(corner_wgt) <- paste('Vol ', round(exp_vol * 100, 1))
round(corner_wgt * 100, 1)
#> Vol 14.5 Vol 14.5 Vol 8.1 Vol 7.4 Vol 7.4 Vol 7.4 Vol 7.4
#> IVV 50 50 9.2 0.7 0.0 0.0 0.0
#> IWM 0 0 0.0 0.0 0.0 0.0 0.0
#> VTV 0 0 0.0 0.0 0.3 0.0 0.0
#> IWF 50 50 50.0 50.0 50.0 50.0 50.0
#> EFA 0 0 0.0 0.0 0.0 0.0 0.0
#> VWO 0 0 0.0 0.0 0.0 0.0 0.0
#> TLT 0 0 40.8 31.1 30.7 30.4 30.3
#> IEF 0 0 0.0 0.0 0.0 0.0 0.0
#> EMB 0 0 0.0 0.0 0.0 0.0 0.0
#> LQD 0 0 0.0 18.2 19.0 19.6 19.7
#> HYG 0 0 0.0 0.0 0.0 0.0 0.0
#> VNQ 0 0 0.0 0.0 0.0 0.0 0.0
#> AMJ 0 0 0.0 0.0 0.0 0.0 0.0
#> IGF 0 0 0.0 0.0 0.0 0.0 0.0
#> GNR 0 0 0.0 0.0 0.0 0.0 0.0
#> GLD 0 0 0.0 0.0 0.0 0.0 0.0
#> Vol 6.8 Vol 6.1 Vol 5 Vol 5 Vol 5 Vol 4.8 Vol 4.8 Vol 4.3
#> IVV 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
#> IWM 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
#> VTV 0.0 0.0 2.0 2.0 2.0 4.4 5.0 6.8
#> IWF 44.3 39.0 28.6 28.6 28.6 25.2 24.2 15.6
#> EFA 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
#> VWO 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
#> TLT 23.9 14.8 0.0 0.0 0.0 0.0 0.0 0.0
#> IEF 0.0 14.8 39.5 39.5 39.5 42.6 43.3 50.0
#> EMB 0.0 0.0 0.0 0.0 0.0 0.0 0.6 2.1
#> LQD 31.8 31.4 29.9 29.9 29.9 27.9 26.9 15.8
#> HYG 0.0 0.0 0.0 0.0 0.0 0.0 0.0 9.7
#> VNQ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
#> AMJ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
#> IGF 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
#> GNR 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
#> GLD 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
#> Vol 4.3 Vol 4.3 Vol 4.3 Vol 4.3 Vol 4 Vol 4
#> IVV 0.0 0.0 0.0 0.0 6.7 0.0
#> IWM 0.0 0.0 0.0 0.0 0.0 0.0
#> VTV 6.8 6.8 6.9 7.4 9.6 13.7
#> IWF 15.6 15.6 15.3 14.3 0.0 0.0
#> EFA 0.0 0.0 0.0 0.0 0.0 0.0
#> VWO 0.0 0.0 0.0 0.0 0.0 0.0
#> TLT 0.0 0.0 0.0 0.0 0.0 0.0
#> IEF 50.0 50.0 50.0 50.0 50.0 50.0
#> EMB 2.1 2.1 2.2 2.8 9.3 12.5
#> LQD 15.8 15.8 15.6 15.5 14.4 13.9
#> HYG 9.7 9.7 10.0 10.0 10.0 10.0
#> VNQ 0.0 0.0 0.0 0.0 0.0 0.0
#> AMJ 0.0 0.0 0.0 0.0 0.0 0.0
#> IGF 0.0 0.0 0.0 0.0 0.0 0.0
#> GNR 0.0 0.0 0.0 0.0 0.0 0.0
#> GLD 0.0 0.0 0.0 0.0 0.0 0.0calc_target_vol solves for the weights for an efficient portfolio with a specific volatility target. If the volatility target is not on the frontier a warning message will print and a value of NA will be returned.
vol_10_wgt <- calc_target_vol(store_constr, cov_mat, 0.1)
rownames(vol_10_wgt) <- asset_class
round(vol_10_wgt * 100, 1)
#> [,1]
#> IVV 24.9
#> IWM 0.0
#> VTV 0.0
#> IWF 50.0
#> EFA 0.0
#> VWO 0.0
#> TLT 25.1
#> IEF 0.0
#> EMB 0.0
#> LQD 0.0
#> HYG 0.0
#> VNQ 0.0
#> AMJ 0.0
#> IGF 0.0
#> GNR 0.0
#> GLD 0.0
sqrt(t(vol_10_wgt) %*% cov_mat %*% vol_10_wgt)
#> [,1]
#> [1,] 0.09999382