Portfolio Optimization with R


Portfolio Optimization with R

According to the Modern portfolio theory (MPT) for any given level of risk it is possible to maximize the return of a portfolio, which is in practice called portfolio optimization. To do this oone would need the historical prices of the assets that will be used to compute mean returns for the time period, as well as the covariance matrix between the assets for the same period, and finally random weights assigned to each asset and to maximize the return to risk ratio.

First you have to install and load the following packages:

library(tidyquant) 
library(timetk) 
library(forcats)
library(tidyr)
library(kableExtra)
library(ggplot2)
library(dplyr)

Sector etfs

This analysis is based on sector ETFs to gain a perspective on the performance and risk of different sectors. IXC energy sector, IXG financial sector, IXN technology sector, IXJ healthcare sector, IXP telecom sector, RXI consumer discretionary sector, EXI industrial sector, MXI basic sector, KXI consumer staple sector, and forJXI utlities sector. Here we use the tq_get function.

tick <- c('IXC', 'IXG', 'IXN', 'IXJ', 'IXP','RXI','EXI','MXI','KXI','JXI')

price_data <- tq_get(tick,
                     from = '2010-01-01',
                     to = '2021-11-01',
                     get = 'stock.prices')

as always we transform the price to return and log transform it using tq_transmute function

log_ret_tidy <- price_data %>%
  dplyr::group_by(symbol) %>%
  tq_transmute(select = adjusted,
               mutate_fun = periodReturn,
               period = 'daily',
               col_rename = 'ret',
               type = 'log')

head(log_ret_tidy)%>%kable()
symbol date ret
IXC 2010-01-04 0.0000000
IXC 2010-01-05 0.0054231
IXC 2010-01-06 0.0086161
IXC 2010-01-07 -0.0029534
IXC 2010-01-08 0.0037575
IXC 2010-01-11 0.0074727

then we transform the long data to wide data using the spread function and drop the missing data using the drop_na function.

log_ret_xts <- log_ret_tidy %>%
  spread(symbol, value = ret) %>%
  tk_xts()
## Warning: Non-numeric columns being dropped: date
## Using column `date` for date_var.
log_ret_xts=log_ret_xts%>%as.data.frame()%>%drop_na()

summary(log_ret_xts)%>%kable()
EXI IXC IXG IXJ IXN IXP JXI KXI MXI RXI
Min. :-0.1079408 Min. :-2.163e-01 Min. :-0.1415391 Min. :-0.0979635 Min. :-0.1537118 Min. :-0.1025307 Min. :-0.1241903 Min. :-0.0967470 Min. :-0.1191473 Min. :-0.1300503
1st Qu.:-0.0046489 1st Qu.:-6.889e-03 1st Qu.:-0.0055047 1st Qu.:-0.0040270 1st Qu.:-0.0048440 1st Qu.:-0.0046377 1st Qu.:-0.0044924 1st Qu.:-0.0034824 1st Qu.:-0.0061805 1st Qu.:-0.0044866
Median : 0.0009071 Median : 4.926e-04 Median : 0.0008703 Median : 0.0008421 Median : 0.0012534 Median : 0.0006725 Median : 0.0007154 Median : 0.0005264 Median : 0.0003481 Median : 0.0009753
Mean : 0.0004047 Mean : 4.899e-05 Mean : 0.0002916 Mean : 0.0004754 Mean : 0.0006594 Mean : 0.0003062 Mean : 0.0002286 Mean : 0.0003485 Mean : 0.0001967 Mean : 0.0005240
3rd Qu.: 0.0062073 3rd Qu.: 7.699e-03 3rd Qu.: 0.0068296 3rd Qu.: 0.0056683 3rd Qu.: 0.0069806 3rd Qu.: 0.0056638 3rd Qu.: 0.0055806 3rd Qu.: 0.0047739 3rd Qu.: 0.0072199 3rd Qu.: 0.0064221
Max. : 0.0990035 Max. : 1.597e-01 Max. : 0.1066680 Max. : 0.0668417 Max. : 0.1050387 Max. : 0.0709797 Max. : 0.1014732 Max. : 0.0786298 Max. : 0.1034600 Max. : 0.0961797
head(log_ret_xts)%>%kable()
EXI IXC IXG IXJ IXN IXP JXI KXI MXI RXI
2010-01-04 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
2010-01-05 0.0017576 0.0054231 0.0123934 -0.0089702 0.0017488 -0.0019752 -0.0024789 -0.0106712 0.0021808 0.0022542
2010-01-06 0.0021924 0.0086161 0.0023331 0.0034450 -0.0071896 -0.0104807 0.0028915 -0.0003517 0.0175849 0.0004501
2010-01-07 0.0074191 -0.0029534 0.0038065 -0.0030617 -0.0026436 -0.0102250 -0.0076597 -0.0035249 -0.0045976 0.0049393
2010-01-08 0.0121002 0.0037575 0.0084070 0.0049704 0.0080843 0.0012839 0.0051818 0.0001765 0.0100872 0.0080304
2010-01-11 0.0123799 0.0074727 -0.0002095 0.0091116 -0.0031557 0.0016483 0.0115111 0.0056329 0.0053080 0.0042126

using the long data colMeans function provides the mean daily return for each sector

mean_ret <- colMeans(log_ret_xts,na.rm = T)
print(round(mean_ret, 5))%>%kable()
##     EXI     IXC     IXG     IXJ     IXN     IXP     JXI     KXI     MXI     RXI 
## 0.00040 0.00005 0.00029 0.00048 0.00066 0.00031 0.00023 0.00035 0.00020 0.00052
x
EXI 0.00040
IXC 0.00005
IXG 0.00029
IXJ 0.00048
IXN 0.00066
IXP 0.00031
JXI 0.00023
KXI 0.00035
MXI 0.00020
RXI 0.00052

in the same way the daily covariance matrix is formed by the cov function and multiplied by the number of work days per year i.e. 252 days excluding weekends and holidays.

cov_mat <- cov(log_ret_xts) * 252
print(round(cov_mat,4))%>%kable()
##        EXI    IXC    IXG    IXJ    IXN    IXP    JXI    KXI    MXI    RXI
## EXI 0.0365 0.0397 0.0390 0.0245 0.0322 0.0254 0.0239 0.0202 0.0379 0.0321
## IXC 0.0397 0.0679 0.0468 0.0267 0.0353 0.0296 0.0269 0.0224 0.0465 0.0355
## IXG 0.0390 0.0468 0.0492 0.0275 0.0358 0.0293 0.0273 0.0224 0.0427 0.0359
## IXJ 0.0245 0.0267 0.0275 0.0251 0.0260 0.0205 0.0197 0.0171 0.0263 0.0237
## IXN 0.0322 0.0353 0.0358 0.0260 0.0418 0.0277 0.0227 0.0207 0.0355 0.0327
## IXP 0.0254 0.0296 0.0293 0.0205 0.0277 0.0271 0.0205 0.0175 0.0285 0.0256
## JXI 0.0239 0.0269 0.0273 0.0197 0.0227 0.0205 0.0282 0.0185 0.0258 0.0222
## KXI 0.0202 0.0224 0.0224 0.0171 0.0207 0.0175 0.0185 0.0193 0.0219 0.0193
## MXI 0.0379 0.0465 0.0427 0.0263 0.0355 0.0285 0.0258 0.0219 0.0481 0.0345
## RXI 0.0321 0.0355 0.0359 0.0237 0.0327 0.0256 0.0222 0.0193 0.0345 0.0339
EXI IXC IXG IXJ IXN IXP JXI KXI MXI RXI
EXI 0.0365 0.0397 0.0390 0.0245 0.0322 0.0254 0.0239 0.0202 0.0379 0.0321
IXC 0.0397 0.0679 0.0468 0.0267 0.0353 0.0296 0.0269 0.0224 0.0465 0.0355
IXG 0.0390 0.0468 0.0492 0.0275 0.0358 0.0293 0.0273 0.0224 0.0427 0.0359
IXJ 0.0245 0.0267 0.0275 0.0251 0.0260 0.0205 0.0197 0.0171 0.0263 0.0237
IXN 0.0322 0.0353 0.0358 0.0260 0.0418 0.0277 0.0227 0.0207 0.0355 0.0327
IXP 0.0254 0.0296 0.0293 0.0205 0.0277 0.0271 0.0205 0.0175 0.0285 0.0256
JXI 0.0239 0.0269 0.0273 0.0197 0.0227 0.0205 0.0282 0.0185 0.0258 0.0222
KXI 0.0202 0.0224 0.0224 0.0171 0.0207 0.0175 0.0185 0.0193 0.0219 0.0193
MXI 0.0379 0.0465 0.0427 0.0263 0.0355 0.0285 0.0258 0.0219 0.0481 0.0345
RXI 0.0321 0.0355 0.0359 0.0237 0.0327 0.0256 0.0222 0.0193 0.0345 0.0339

The runif function provides hypothetical weights standardized to sum up to 1.

wts <- runif(n = length(tick))
wts <- wts/sum(wts)
print(wts)%>%kable()
##  [1] 0.02181918 0.07788497 0.05416218 0.21472448 0.06099253 0.04694185
##  [7] 0.19240680 0.15534889 0.07119757 0.10452155
x
0.0218192
0.0778850
0.0541622
0.2147245
0.0609925
0.0469418
0.1924068
0.1553489
0.0711976
0.1045216

This is the code for computing portfolio returns, risks, and sharp ratios.

port_returns <- (sum(wts * mean_ret) + 1)^252 - 1


port_risk <- sqrt(t(wts) %*%(cov_mat %*% wts) )
print(port_risk)%>%kable()
##          [,1]
## [1,] 0.158721
0.158721
sharpe_ratio <- port_returns/port_risk
print(sharpe_ratio)%>%kable()
##           [,1]
## [1,] 0.5842914
0.5842914

This snippet set place holders for 10000 weights, portfolio risks, returns, and sharp ratios.

num_port <- 10000

# Creating a matrix to store the weights

all_wts <- matrix(nrow = num_port,
                  ncol = length(tick))

# Creating an empty vector to store
# Portfolio returns

port_returns <- vector('numeric', length = num_port)

# Creating an empty vector to store
# Portfolio Standard deviation

port_risk <- vector('numeric', length = num_port)

# Creating an empty vector to store
# Portfolio Sharpe Ratio

sharpe_ratio <- vector('numeric', length = num_port)

here is the loop to highlight portfolio returns, risks, and sharp ratios for different wheight combinations

for (i in seq_along(port_returns)) {
  
  wts <- runif(length(tick))
  wts <- wts/sum(wts)
  
  # Storing weight in the matrix
  all_wts[i,] <- wts
  
  # Portfolio returns
  
  port_ret <- sum(wts * mean_ret)
  port_ret <- ((port_ret + 1)^252) - 1
  
  # Storing Portfolio Returns values
  port_returns[i] <- port_ret
  
  
  # Creating and storing portfolio risk
  port_sd <- sqrt(t(wts) %*% (cov_mat  %*% wts))
  port_risk[i] <- port_sd
  
  # Creating and storing Portfolio Sharpe Ratios
  # Assuming 0% Risk free rate
  
  sr <- port_ret/port_sd
  sharpe_ratio[i] <- sr
  
}

here we are making a tibble of portfolio returns, risks, and sharp ratios

portfolio_values <- tibble(Return = port_returns,
                           Risk = port_risk,
                           SharpeRatio = sharpe_ratio)


# Converting matrix to a tibble and changing column names
all_wts <- tk_tbl(all_wts)
## Warning in tk_tbl.data.frame(as.data.frame(data), preserve_index,
## rename_index, : Warning: No index to preserve. Object otherwise converted to
## tibble successfully.
colnames(all_wts) <- colnames(log_ret_xts)

# Combing all the values together
portfolio_values <- tk_tbl(cbind(all_wts, portfolio_values))
## Warning in tk_tbl.data.frame(cbind(all_wts, portfolio_values)): Warning: No
## index to preserve. Object otherwise converted to tibble successfully.
head(portfolio_values)%>%kable()
EXI IXC IXG IXJ IXN IXP JXI KXI MXI RXI Return Risk SharpeRatio
0.0015059 0.1415284 0.1188994 0.1288394 0.0145800 0.1913636 0.1803470 0.0978158 0.0370134 0.0881072 0.0786358 0.1649985 0.4765849
0.1486655 0.0318730 0.0683066 0.0787044 0.1056284 0.0336204 0.1851120 0.1775181 0.0860809 0.0844906 0.0961518 0.1626657 0.5911005
0.0671929 0.1721032 0.0626337 0.0954884 0.1794551 0.0124373 0.0973874 0.1729593 0.0886350 0.0517077 0.0917700 0.1720987 0.5332403
0.1399838 0.0752763 0.0430130 0.0601343 0.0910031 0.1037290 0.1257542 0.0925847 0.1396637 0.1288578 0.0919778 0.1701907 0.5404394
0.1707087 0.0214482 0.0324814 0.1145488 0.1880701 0.0873300 0.0122085 0.0890225 0.1952464 0.0889353 0.1070422 0.1738556 0.6156959
0.1134023 0.0345157 0.1007384 0.0136946 0.1229113 0.1539289 0.0585701 0.1751657 0.1650301 0.0620429 0.0924751 0.1694699 0.5456725

Minimum variance portfolio weights and Tangency portfolio weights with the highest sharp ratio!

min_var <- portfolio_values[which.min(portfolio_values$Risk),]
max_sr <- portfolio_values[which.max(portfolio_values$SharpeRatio),]

Minimum variance portfolio weights visualization is presented here

p <- min_var %>%
  gather(EXI:RXI, key = Asset,
         value = Weights) %>%
  mutate(Asset = as.factor(Asset)) %>%
  ggplot(aes(x = fct_reorder(Asset,Weights), y = Weights, fill = Asset)) +
  geom_bar(stat = 'identity') +
  theme_minimal() +
  labs(x = 'Assets', y = 'Weights', title = "Minimum Variance Portfolio Weights") +
  scale_y_continuous(labels = scales::percent) 

p

Tangency portfolio weights with the highest sharp ratio visualization is presented here

p <- max_sr %>%
  gather(EXI:RXI, key = Asset,
         value = Weights) %>%
  mutate(Asset = as.factor(Asset)) %>%
  ggplot(aes(x = fct_reorder(Asset,Weights), y = Weights, fill = Asset)) +
  geom_bar(stat = 'identity') +
  theme_minimal() +
  labs(x = 'Assets', y = 'Weights', title = "Tangency Portfolio Weights") +
  scale_y_continuous(labels = scales::percent) 

p

and the risk/performance axis for all weights is presented here

p <- portfolio_values %>%
  ggplot(aes(x = Risk, y = Return, color = SharpeRatio)) +
  geom_point() +
  theme_classic() +
  scale_y_continuous(labels = scales::percent) +
  scale_x_continuous(labels = scales::percent) +
  labs(x = 'Annualized Risk',
       y = 'Annualized Returns',
       title = "Portfolio Optimization & Efficient Frontier") +
  geom_point(aes(x = Risk,
                 y = Return), data = min_var, color = 'red') +
  geom_point(aes(x = Risk,
                 y = Return), data = max_sr, color = 'red')
p

References