Quantitative Finance: Bayesian Portfolio Selection

R Finance Series

Challange of Classical Portfolio Selection

Markowitz’s 1952 paper set the foundations for what is now popularly referred to as Modern Portfolio Theory (MPT) and had a profound impact on the financial industry. Individual security selection lay at the heart of the standard investment practice until then. Afterward, the focus shifted toward diversification and assessment of the contribution of individual securities to the risk-return profile of a portfolio. Mean-variance analysis rests on the assumption that rational investors select among risky assets on the basis of expected portfolio return and risk (as measured by the portfolio variance).

The classical approach to asset allocation is a two-step process: In the first step one estimates the parameters for return distribution. In the second stage, portfolio weights using the estimated parameters can be optimized. However, the reputation of this classical approach among practitioners has suffered due to numerous implementation difficulties. Portfolio weights derived from it are notoriously sensitive to the inputs, especially expected returns, and often represent unintuitive or extreme allocations exposing an investor to unintended risks. These inputs (expected returns, variances, and covariances) are all subject to estimation errors that an optimizer picks up and then leverages.

Bayesian Approach for Portfolio Allocation

Bayesian approach provides a way to limit the sensitivity of the final allocation to the input parameters by shrinking the estimate of the market parameters toward the investor’s prior. Using the Bayesian approach to estimation we can naturally identify a suitable uncertainty range for the market parameters, namely the location-dispersion ellipsoid of their posterior distribution. Bayesian allocations are the solutions to a robust optimization problem that uses as uncertainty range the Bayesian location-dispersion ellipsoid.

R Codes and Empirical Results

#===================================================
#    Bayesian Approach for Portfolio Selection
#===================================================

#--------------------
#    Prepare data
#--------------------

rm(list = ls())
library(tidyverse)
library(tidyquant)

symbols <- c("FB", "AMZN", "NFLX", "GOOG", "^IXIC")

tq_get(x = symbols, from = "2017-01-01", to = "2020-12-31") -> data_raw

data_raw %>% 
  select(date, symbol, adjusted) %>% 
  pivot_wider(names_from = symbol, values_from = adjusted) %>% 
  slice(-which.min(date)) -> data_wide_price

date_ymd <- data_wide_price$date
market_data <- data_wide_price$`^IXIC`
stocks_data <- data_wide_price %>% select(symbols[1:4])

data_raw %>%
  filter(symbol %in% symbols[1:4]) %>% 
  group_by(symbol) %>%
  tq_transmute(select = adjusted, 
               mutate_fun = periodReturn, 
               period = "daily", 
               type = "log",
               col_rename = "daily_return") %>% 
  pivot_wider(names_from = symbol, values_from = daily_return) %>% 
  slice(-which.min(date)) -> data_wide_return

#--------------------------------------------------------------
#  Bayesian Approach for Portfolio Selection/Asset Allocation
#--------------------------------------------------------------

# Remove date column: 

data_wide_return %>% select(-date) -> return_wide_train

# Calculate stock mean: 
mean_ret <- colMeans(return_wide_train) 

# Covariance matrix: 

cov_mat <- cov(return_wide_train)

# Number of observations: 
T_rows <- nrow(return_wide_train)

# Number of stock symbols: 
N <- ncol(return_wide_train)

# Set parameters as proposed by Rachev et al. (2008) at page 106, Baysian Method in Finance: 

nu <- N + 2 # (v parameter)
tau <- 500 # (tau parameter)

data_stan <- list(
  T_rows = T_rows,
  N = N,
  nu = nu,
  tau = tau,
  eta = mean_ret,
  R = as.matrix(return_wide_train),
  omega = cov_mat * (nu - N -1)
)

# Stan codes (download from https://www.mediafire.com/file/7q7czce8dtrx56j/my_stan_codes_bayesian_por.stan/file): 
data {
  //5 is arbitrary
  int<lower = 5> T_rows;
  int<lower = 2> N;
  real<lower = N - 1> nu;
  real<lower = 0> tau;
  vector[N] eta;
  matrix[T_rows, N] R;
  cov_matrix[N] omega;
}

parameters {
  vector[N] mu;
  cov_matrix[N] sigma;
}

transformed parameters{
  cov_matrix[N] sigma_scaled;
  sigma_scaled = (1 / tau) * sigma;
}

model {
  
  target += inv_wishart_lpdf(sigma | nu, omega);
  target += multi_normal_lpdf(mu | eta, sigma_scaled);
  for(t in 1:T_rows){
    target += multi_normal_lpdf(R[t]| mu, sigma);
  }
  
}

# Fitting the model: 

library(rstan)

options(mc.cores = parallel::detectCores())

rstan_options(auto_write = TRUE)

# Set the number of Markov chains: 
chains <- N + 2

# Set the number of iterations for each chain: 
iter <- 2000

# No-U-Turn sampler variant of Hamiltonian Monte Carlo 
# (Hoffman and Gelman 2011, Betancourt 2017) process: 

my_algorithm <- "NUTS"

# Fit Bayesian Process: 

system.time(
  fit <- stan(
    file = "C:/Users/ADMIN/Documents/my_stan_codes_bayesian_por.stan",
    data = data_stan,
    algorithm = my_algorithm, 
    chains = chains,
    warmup = iter / 2,
    iter = iter,
    seed = 29, 
    verbose = TRUE
  )
)


list_of_draws <- rstan::extract(fit)

# saveRDS(list_of_draws, file = "list_of_draws.Rds")

# Extract the posterior distribution of sigma: 

sigma_post <- list_of_draws$sigma

# Extract the posterior distribution of mu: 

mu_post <- list_of_draws$mu

# Function extracts weights by Bayesian Approach: 

weights_bayesian <- function(i) {
  
  bayesian_sigma_port <- sigma_post[i, , ]
  
  top_mat <- cbind(2*bayesian_sigma_port, rep(1, N))
  
  bot_vec <- c(rep(1, N), 0)
  
  Am_mat <- rbind(top_mat, bot_vec)
  
  b_vec <- c(rep(0, N), 1)
  
  zm_mat <- solve(Am_mat) %*% b_vec
  
  X_min_var_baysian <- zm_mat[1:N, 1]
  
  return(X_min_var_baysian)
}


# Calculate Bayesian weights: 

number_sigmas <- length(sigma_post) / (N*N)

sapply(1:number_sigmas, weights_bayesian) %>% 
  t() %>% 
  colMeans() -> bayesian_weights

# Bayesian approach for asset allocation: 

bayesian_portfolio <- function(...) {
  
  bayesian_sigma_port <- sigma_post[1, , ]
  
  top_mat <- cbind(2*bayesian_sigma_port, rep(1, N))
  
  bot_vec <- c(rep(1, N), 0)
  
  Am_mat <- rbind(top_mat, bot_vec)
  
  b_vec <- c(rep(0, N), 1)
  
  zm_mat <- solve(Am_mat) %*% b_vec
  
  X_min_var_baysian <- zm_mat[1:N, 1]
  
  return(X_min_var_baysian)
}


# Traditional approach (as proposed by Markowitz) for asset allocation: 

markowitz_portfolio <- function(...) {
  
  top_mat <- cbind(2*cov_mat, rep(1, N))
  
  bot_vec <- c(rep(1, N), 0)
  
  Am_mat <- rbind(top_mat, bot_vec)
  
  b_vec <- c(rep(0, N), 1)
  
  zm_mat <- solve(Am_mat) %*% b_vec
  
  X_min_var_markowitz <- zm_mat[1:N, 1]
  
  return(X_min_var_markowitz)

}


#------------------------------------------------------------
#  Compare the two methods by Portfolio Backtesting Process
#------------------------------------------------------------


portfolios <- list("Bayesian" = bayesian_portfolio, "Markowitz" = markowitz_portfolio)

library(xts)

xts(stocks_data, order.by = date_ymd) -> stocks_data
xts(market_data, order.by = date_ymd) -> market_data


data <- list(list(adjusted = stocks_data, index = market_data))

library(portfolioBacktest)
bt <- portfolioBacktest(portfolios, data, benchmark = c("uniform", "index"), show_progress_bar = TRUE)

res_sum <- backtestSummary(bt)

res_sum$performance_summary[1:9, 1:2] %>% 
  as.data.frame() %>% 
  mutate(Criterion = row.names(.)) %>% 
  mutate_if(is.numeric, function(x) {round(x, 3)}) %>% 
  select(Criterion, everything()) -> df_compare


library(kableExtra) # For presenting table. 

df_compare %>% 
  kbl(caption = "Table 1: Portfolio Performance by Approach", escape = TRUE) %>%
  kable_classic(full_width = FALSE, html_font = "Cambria")

Table 1: Portfolio Performance by Approach

Criterion	Bayesian	Markowitz
Sharpe ratio	0.915	0.877
max drawdown	0.265	0.269
annual return	0.267	0.256
annual volatility	0.292	0.292
Sterling ratio	1.009	0.951
Omega ratio	1.189	1.183
ROT (bps)	1948.218	1900.753
VaR (0.95)	0.031	0.031
CVaR (0.95)	0.045	0.045

# backtestChartCumReturns(bt, c("Bayesian", "Markowitz"))

bt$Bayesian$data1$return %>% cumsum() %>% 
  as.data.frame() %>% 
  mutate(date_ymd = row.names(.) %>% ymd()) %>% 
  rename(CumReturn = `portfolio return`) %>% 
  mutate(Approach = "Bayesian") -> df1

bt$Markowitz$data1$return %>% cumsum() %>% 
  as.data.frame() %>% 
  mutate(date_ymd = row.names(.) %>% ymd()) %>% 
  rename(CumReturn = `portfolio return`) %>% 
  mutate(Approach  = "Markowitz") -> df2

bind_rows(df1, df2) -> df

df %>% 
  mutate(CumReturn = round(CumReturn, 3), Date = date_ymd) %>% 
  ggplot(aes(Date, CumReturn, color = Approach)) + 
  geom_line() + 
  labs(title = "Figure 1: Daily Cummulative Return by Approach") -> p

plotly::ggplotly(p)