Downloading data

First lets load our packages

library(quantmod)
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(tidyquant) # To download the data
## Loading required package: lubridate
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
## Loading required package: PerformanceAnalytics
## 
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
## ══ Need to Learn tidyquant? ════════════════════════════════════════════════════
## Business Science offers a 1-hour course - Learning Lab #9: Performance Analysis & Portfolio Optimization with tidyquant!
## </> Learn more at: https://university.business-science.io/p/learning-labs-pro </>
library(plotly) # To create interactive charts
## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(timetk) # To manipulate the data series
library(tidyr)
library(ggplot2)
library(forcats)

Next lets select a few stocks to build our portfolios.

We will choose the following 4 stocks

Lets download the price data.

tick <- c('AMZN', 'AAPL', 'NFLX', 'T')

price_data = tq_get(tick,
                    from = '2014-01-01',
                    to = '2022-03-31',
                    get = 'stock.prices',
                    complete_cases=T)

Next we will calculate the daily returns for these stocks. We will use the logarithmic returns.

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

As you can see that this data is in tidy format. We will use the spread() function to convert it to a wide format. And we will also convert it into a time series object using xts() 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.

Next let’s calculate the mean daily returns for each asset. Calculate the covariance matrix for all these stocks. We will annualize it by multiplying by 252.

mean_ret <- colMeans(log_ret_xts)
print(round(mean_ret, 4))
##   AAPL   AMZN   NFLX      T 
## 0.0011 0.0010 0.0010 0.0001
cov_mat <- cov(log_ret_xts) * 252
print(round(cov_mat,4))
##        AAPL   AMZN   NFLX      T
## AAPL 0.0799 0.0456 0.0466 0.0191
## AMZN 0.0456 0.0959 0.0663 0.0126
## NFLX 0.0466 0.0663 0.1787 0.0131
## T    0.0191 0.0126 0.0131 0.0447

To calculate the portfolio returns and risk (standard deviation) we will us need

# Calculate the random weights
wts <- runif(n = length(tick))
wts <- wts/sum(wts)
print(wts)
## [1] 0.1095581 0.2065660 0.2932773 0.3905987
# Calculate the portfolio returns
port_returns <- (sum(wts * mean_ret) + 1)^252 - 1

# Calculate the portfolio risk
port_risk <- sqrt(t(wts) %*% (cov_mat %*% wts))

# Calculate the Sharpe Ratio
sharpe_ratio <- port_returns/port_risk

We have everything we need to perform our optimization. All we need now is to run this code on 5000 random portfolios. For that we will use a for loop.

Before we do that, we need to create empty vectors and matrix for storing our values.

# Creating a matrix to store the weights
num_port <- 10000
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)

Next lets run the for loop 5000 times.

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
}

All the heavy lifting has been done and now we can create a data table to store all the values together.

# Storing the values in the table
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.

We have the weights in each asset with the risk and returns along with the Sharpe ratio of each portfolio.

Next lets look at the portfolios that matter the most.

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

Plot the weights of each portfolio. First with the minimum variance portfolio.

p <- min_var %>%
  gather(AAPL:T, 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) 

ggplotly(p)

Next lets look at the tangency portfolio or the the portfolio with the highest sharpe ratio.

p <- max_sr %>%
  gather(AAPL:T, 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) 

ggplotly(p)

Finally lets plot all the random portfolios and visualize the efficient frontier.

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') +
  annotate('text', x = 0.20, y = 0.42, label = "Tangency Portfolio") +
  annotate('text', x = 0.18, y = 0.01, label = "Minimum variance portfolio") +
  annotate(geom = 'segment', x = 0.14, xend = 0.135,  y = 0.01, 
           yend = 0.06, color = 'red', arrow = arrow(type = "open")) +
  annotate(geom = 'segment', x = 0.22, xend = 0.2275,  y = 0.405, 
           yend = 0.365, color = 'red', arrow = arrow(type = "open"))
  

ggplotly(p)