portfolio_optimization.md

title	author
Portfolio Optimization & Rebalancing Tool	John Houghton

Setup R Environment

library(knitr)
library(rmdformats)
## Global options
options(max.print="75")
opts_chunk$set(#echo=FALSE,
               cache=TRUE,
               prompt=FALSE,
               tidy=TRUE,
               comment=NA,
               message=FALSE,
               warning=FALSE, 
               results='asis'
               )
opts_knit$set(width=75)
print("notebook options imported successfully!")

## [1] "notebook options imported successfully!"

## Load Required Libraries & Set Paths
library(tidyverse)
library(tidyquant)
library(lubridate)
library(timetk)
library(DT)
library(PerformanceAnalytics)
library(ROI.plugin.quadprog)
library(ROI.plugin.glpk)
# library(ROI.plugin.symphony)
setwd("~/Documents/GitHub/rob_gordon_2018")

Get Stock Data

stocks <- c("SPY", "IYR", "QQQ", "TLT", "GLD", "IWM", "EEM", "DBC", "EFA")
stock_data <- tq_get(stocks, get = "stock.prices", from = Sys.Date() - months(12), 
    to = Sys.Date())
stock_data %>% head(10)

A tibble: 10 x 8

symbol date open high low close volume adjusted 1 SPY 2017-11-20 258. 259. 258. 258. 48075500 254. 2 SPY 2017-11-21 259. 260. 258. 260. 69176800 255. 3 SPY 2017-11-22 260 260. 260. 260. 45033400 255. 4 SPY 2017-11-24 260. 260. 260. 260. 27856500 256. 5 SPY 2017-11-27 260. 261. 260 260. 52274900 256. 6 SPY 2017-11-28 261. 263. 261. 263. 98971700 258. 7 SPY 2017-11-29 263. 264. 262. 263. 77512100 258. 8 SPY 2017-11-30 264. 266. 264. 265. 127894400 260. 9 SPY 2017-12-01 265. 265. 261. 264. 164390900 260. 10 SPY 2017-12-04 266. 267. 264. 264. 94040600 259.

Visualize Stock Prices

• Kept only 4 symbols in order to fit
• Will include all in dashboard via interactive filtering

# using ggplot2
stock_data %>% filter(symbol %in% c("SPY", "GLD", "QQQ", "DBC")) %>% ggplot(aes(x = date, 
    y = adjusted, color = symbol)) + geom_line(size = 1) + geom_bbands(aes(high = high, 
    low = low, close = close), ma_fun = SMA, n = 5, show.legend = TRUE) + labs(title = "Daily Stock Prices", 
    x = "", y = "Adjusted Prices", color = "") + facet_wrap(~symbol, ncol = 2, 
    scales = "free") + scale_y_continuous(labels = scales::dollar) + theme_tq() + 
    scale_color_tq()

Calculate Returns{.tabset .tabset-fade}

Simple Returns Data

mo_returns <- stock_data %>% group_by(symbol) %>% tq_transmute(select = adjusted, 
    mutate_fun = periodReturn, period = "monthly", col_rename = "returns")

mo_returns %>% mutate(returns = round(returns, 6)) %>% arrange(desc(symbol), 
    desc(date)) %>% datatable(class = "compact", options = list(scrollX = TRUE, 
    dom = "tp"))

Error in loadNamespace(name): there is no package called 'webshot'

# simple returns = Rt = [ Pt - Pt-1 ]/ Pt-1 simple_returns <- stock_data %>%
# group_by(symbol, year = year(date)) %>% summarise(adjusted =
# last(adjusted, order_by = year)) %>% mutate(return = (adjusted -
# lag(adjusted, order_by = year)) / lag(adjusted, order_by = year) , return
# = round(return, 6)) %>% arrange(desc(symbol), desc(year)) # %>%
# filter(symbol == 'SPY') datatable(simple_returns, class='compact', options
# = list(scrollX = TRUE, dom = 'tp'))

Chart: Returns

• stocks filtered to fit

mo_returns %>% filter(symbol %in% c("SPY", "GLD", "QQQ", "DBC")) %>% ggplot(aes(x = date, 
    y = returns, fill = symbol)) + geom_line(aes(color = symbol)) + geom_hline(yintercept = 0, 
    color = palette_light()[[1]]) + scale_y_continuous(labels = scales::percent) + 
    labs(title = " Monthly Returns", subtitle = "stocks limited to fit", y = "Monthly Returns", 
        x = "") + facet_wrap(~symbol, ncol = 2, scales = "free") + theme_tq() + 
    scale_fill_tq()

Chart: Smoothed Returns

mo_returns %>% filter(symbol %in% c("SPY", "GLD", "QQQ", "DBC")) %>% ggplot(aes(x = date, 
    y = returns, fill = symbol)) + geom_smooth(aes(color = symbol), method = "loess", 
    formula = y ~ x, se = FALSE) + geom_hline(yintercept = 0, color = palette_light()[[1]]) + 
    scale_y_continuous(labels = scales::percent) + labs(title = " Monthly Returns", 
    subtitle = "stocks limited to fit", y = "Monthly Returns", x = "") + facet_wrap(~symbol, 
    ncol = 2, scales = "free") + theme_tq() + scale_fill_tq()

Growth of $1000 (6mo){.tabset .tabset-fade}

Data Table

init.investment <- 1000
growth <- mo_returns %>% arrange(date) %>% mutate(final_value = init.investment * 
    cumprod(1 + returns)) %>% arrange(desc(final_value))
growth %>% filter(date == max(date)) %>% select(-date)

A tibble: 9 x 3

Groups: symbol [9]

symbol returns final_value 1 QQQ -0.0457 1062. 2 SPY -0.00565 1061. 3 IYR 0.0325 1022. 4 IWM -0.00780 1008. 5 DBC -0.0495 996. 6 GLD 0.00452 954. 7 TLT 0.0151 935. 8 EFA -0.000800 931. 9 EEM 0.0304 876.

Visualization

growth %>% ggplot(aes(x = date, y = final_value, color = symbol)) + geom_line() + 
    # geom_smooth(method = 'loess') +
labs(title = "Individual Portfolio: Comparing the Growth of $1K", subtitle = "Quickly visualize performance", 
    x = "", y = "Investment Value") + theme_tq() + theme(legend.position = "right") + 
    scale_y_continuous(labels = scales::dollar)

winners (top 5)

winners <- growth %>% ungroup() %>% filter(date == max(date)) %>% mutate(rank = row_number()) %>% 
    top_n(5, final_value) %>% select(rank, symbol, final_value)

winners %>% arrange(desc(final_value)) %>% ggplot(aes(x = symbol, y = final_value, 
    fill = symbol)) + geom_bar(stat = "identity") + # geom_smooth(method = 'loess') +
labs(title = "Individual Portfolio: Comparing the Growth of $1K", subtitle = "Quickly visualize performance", 
    x = "", y = "Investment Value") + theme_tq() + theme(legend.position = "right") + 
    scale_y_continuous(labels = scales::dollar) + geom_text(aes(label = round(final_value, 
    2), vjust = 3, size = 10))

winners$symbol

[1] "QQQ" "SPY" "IYR" "IWM" "DBC"

Portfolio Optimization

Create Time-Series Object for Optimization

# convert to ts object
stock_returns <- mo_returns %>% filter(symbol %in% winners$symbol) %>% spread(key = symbol, 
    value = returns) %>% tk_xts(date_var = date) %>% na.omit()

# stock_returns %>% class()

Optimize & BACKTEST

Get Historical Stock Returns (3 Years)

hist_returns <- tq_get(winners$symbol, get = "stock.prices", from = Sys.Date() - 
    months(36), to = Sys.Date()) %>% group_by(symbol) %>% tq_transmute(select = adjusted, 
    mutate_fun = periodReturn, period = "daily", type = "arithmetic", col_rename = "returns")

hist_returns_ts <- hist_returns %>% spread(key = symbol, value = returns) %>% 
    tk_xts(date_var = date) %>% na.omit()

Run Minimum Variance Optimization

• set objective: risk minimization (variance)
• contraint 1: invest all funds
• contrstaint 2: long only

# Markowitz Optimization involves minimizing the risk(variance) of returns and maximizing returns.
# https://rdrr.io/cran/PortfolioAnalytics/f/inst/doc/ROI_vignette.pdf
library(PortfolioAnalytics)

# Loading all suggested packages to fix ROI issue
# list.of.packages <- c("foreach", "DEoptim", "iterators", "fGarch", "Rglpk", "quadprog", "ROI", "ROI.plugin.glpk", "ROI.plugin.quadprog", "ROI.plugin.symphony", "pso", "GenSA", "corpcor", "testthat", "nloptr", "MASS", "robustbase")
# new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
# if(length(new.packages)) install.packages(new.packages)

#' Set up initial portfolio with basic constraints.
init.portf <- portfolio.spec(assets = colnames(hist_returns_ts))
# Add full investment constraint to the portfolio object
init.portf <- add.constraint(portfolio=init.portf, type="full_investment")
# Add long only constraints
init.portf <- add.constraint(portfolio=init.portf, type="box", min_sum=0.99, max_sum=1.01)

# Add objective to minimize variance
portf_minvar <- add.objective(portfolio=init.portf, type="risk", name="var")
# Dummy objectives for plotting and/or further analysis
portf_minvar <- add.objective(portf_minvar, type="return", name="mean", multiplier=0)
portf_minvar <- add.objective(portf_minvar, type="risk", name="ES", multiplier=0)
portf_minvar <- add.objective(portf_minvar, type="risk", name="StdDev", multiplier=0)

opt_gmv <- optimize.portfolio(R=hist_returns_ts,
                              portfolio=portf_minvar,
                              optimize_method="ROI",
                              trace=TRUE,
                              search_size=5000)

### Print Min Variance Achieved via Optimization ###
extractStats(opt_gmv)[1]

    ES 

0.01513927

### Print Weights ###
weights <- extractWeights(opt_gmv)
weights_tbl <- as.tibble(as.list(weights)) %>% gather("symbol", "weight", 1:length(weights)) %>% arrange(desc(weight))
weights_tbl

A tibble: 5 x 2

symbol weight 1 IYR 0.420 2 DBC 0.366 3 SPY 0.214 4 IWM 0
5 QQQ 0

### Bar Chart Weights ###
opt_gmv %>% chart.Weights(plot.type="barplot", ylim=c(0,1), las=1, colorset = rainbow12equal, main="Optimal Weights", cex.axis = .8 )

### Pie Chart Weights ###
library(scales)
pie <- ggplot(weights_tbl, aes(x="", y=weight, fill=symbol)) +
  geom_bar(width = 1, stat = "identity") + 
  coord_polar("y", start=0) +
  scale_fill_brewer("Blues") + #blank_theme() +
  theme(axis.text.x=element_blank())+
  geom_text(aes(y = weight/3 + c(0, cumsum(weight)[-length(weight)]), 
                label = percent(weight)), size=5)
pie

paste("Min Variance Achieved: ", round(opt_gmv$objective_measures$ES,5))

[1] "Min Variance Achieved: 0.01514"

### Chart Risk Reward ### 
chart.RiskReward(opt_gmv, rp=TRUE, risk.col = "ES", chart.assets=TRUE)

annualized.moments <- function(R, scale=12, portfolio=NULL){
    out <- list()
    out$mu <-    matrix(colMeans(R), ncol=1)
    out$sigma <- cov(R)/scale
    return(out)
  }

### CALCULATE & PLOT EFFICIENT FRONTIER ### 
  prt_ef <- create.EfficientFrontier(R=12*hist_returns_ts, portfolio=portf_minvar, type="mean-StdDev", 
                                      match.col = "StdDev", momentFUN="annualized.moments", scale=12)

Error in get(as.character(FUN), mode = "function", envir = envir): object 'annualized.moments' of mode 'function' was not found

  ### Chart Efficient Frontier ### 
  xlim <- range(prt_ef$frontier[,2])*c(1, 1.5)

Error in eval(expr, envir, enclos): object 'prt_ef' not found

  ylim <- range(prt_ef$frontier[,1])*c(.80, 1.05)

Error in eval(expr, envir, enclos): object 'prt_ef' not found

  chart.EfficientFrontier(prt_ef, match.col="StdDev", chart.assets = TRUE, n.portfolios = 1000
                          , labels.assets = TRUE, xlim=xlim, ylim=ylim, main = "Efficient Frontier")

Error in chart.EfficientFrontier(prt_ef, match.col = "StdDev", chart.assets = TRUE, : object 'prt_ef' not found

  points(with(annualized.moments(12*stock_returns, scale=12), cbind(sqrt(diag(sigma)), mu)), pch=19 ) 
  text(with(annualized.moments(12*stock_returns, scale=12), cbind(sqrt(diag(sigma)), mu)), 
       labels=colnames(stock_returns), cex=.8, pos=4) 

  ### Chart EF Weights ###
  chart.EF.Weights(prt_ef, match.col="StdDev", colorset = rainbow12equal, main = "Efficient Frontier")

Error in chart.EF.Weights(prt_ef, match.col = "StdDev", colorset = rainbow12equal, : object 'prt_ef' not found

Visualize Portfolio Optimization

Perform Backtest

# https://rdrr.io/cran/PortfolioAnalytics/f/inst/doc/ROI_vignette.pdf

bt_gmv <- optimize.portfolio.rebalancing(R=hist_returns_ts, portfolio=portf_minvar,
                                         optimize_method="ROI",
                                         trace = TRUE,
                                         rebalance_on="months",
                                         training_period=36,
                                         # trailing_period=trailing, 
                                         rolling_window = 12
                                         )

summary(bt_gmv)$annualized_returns[[1]]

[1] 0.1133622

summary(bt_gmv)$annualized_StdDev[[1]]

[1] 0.1018277

bt_weights <- summary(bt_gmv)$weights
bt_ret <- summary(bt_gmv)$portfolio_returns

charts.PerformanceSummary(cbind(bt_ret, opt_gmv$R), main="Optimization Performance", method = "StdDev", wealth.index = TRUE, begin = "axis")

x <- Return.portfolio(hist_returns_ts, bt_weights, wealth.index = TRUE, verbose = TRUE)
chart.CumReturns(x$returns, col=bluemono, begin="axis")

chart.StackedBar(x$BOP.Weight)

chart.StackedBar(x$BOP.Value)

# CHART BACKTEST WEIGHTS
bt_gmv %>% chart.Weights(ylim=c(0,1), las=1, main="Optimal Weights", cex.axis = 3 ,cex.legend = 3, col=rainbow12equal) # , legend.loc = "topright"

# ret.bt.gmv <- do.call(cbind, lapply(bt_gmv, function(x) summary(bt_gmv)$portfolio_returns))
# colnames(ret.bt.gmv) <- c("min ES", "min ES RB", "min ES Eq RB", "stuff", "other stuff")

Rb <- Return.portfolio(hist_returns_ts, weights = c(.2,.2,.2,.2,.2), wealth.index = TRUE, verbose = TRUE)
chart.RelativePerformance(as.xts(x$returns), as.xts(Rb$returns), colorset = rainbow12equal, main = "Performance Relative to Even Dist")

# `r library(dplyr, PerformanceAnalytics)` # Get Performance {.tabset
# .tabset-fade} ## descriptive statistics `r mo_returns %>%
# tidyquant::tq_performance(Ra = returns, Rb = NULL, performance_fun =
# table.Stats)` ## annualized returns `r mo_returns %>%
# tidyquant::tq_performance(Ra = returns, Rb = NULL, performance_fun =
# table.AnnualizedReturns)` ## downside risks `r mo_returns %>%
# tidyquant::tq_performance(Ra = returns, Rb = NULL, performance_fun =
# table.DownsideRisk)` ## downside risk ratios `r mo_returns %>%
# tidyquant::tq_performance(Ra = returns, Rb = NULL, performance_fun =
# table.DownsideRiskRatio)` ## value at risk `r mo_returns %>%
# tidyquant::tq_performance(Ra = returns, Rb = NULL, performance_fun = VaR)`
# ## sharpe ratio `r mo_returns %>% tidyquant::tq_performance(Ra = returns,
# Rb = NULL, performance_fun = SharpeRatio)`

Appendix

Resources & Notes

• Performance Analysis Vignette
• Really great article I walked through
• Really great article #2 - same as above or extended?
• Backtesting 4 Portfolio Optimization Strategies In R
• Using ROI Solvers for Portfolio Analytics, Vignette
• Example R Code from Rob via Systematic Investor
• Systematic Investor Github with Library
• Article from Rob re: strategy
• Moments Lookback WIndow Optimization

Second EF Example from Stack Overflow

require(quadprog)
# min_x(-d^T x + 1/2 b^T D x) r.t A.x>=b
MV_QP <- function(nx, tarRet, Sig = NULL, long_only = FALSE) {
    if (is.null(Sig)) 
        Sig = cov(nx)
    dvec = rep(0, ncol(Sig))
    meq = 2
    Amat = rbind(rep(1, ncol(Sig)), apply(nx, 2, mean))
    bvec = c(1, tarRet)
    if (long_only) {
        meq = 1
        Amat = Amat[-1, ]
        Amat = rbind(Amat, diag(1, ncol(Sig)), rep(1, ncol(Sig)), rep(-1, ncol(Sig)))
        bvec = bvec[-1]
        bvec = c(bvec, rep(0, ncol(Sig)), 0.98, -1.02)
    }
    sol <- solve.QP(Dmat = Sig, dvec, t(Amat), bvec, meq = meq)$solution
}

steps = 50
x = stock_returns
µ.b <- apply(X = x, 2, FUN = mean)
long_only = TRUE
range.bl <- seq(from = min(µ.b), to = max(µ.b) * ifelse(long_only, 1, 1.6), 
    length.out = steps)
risk.bl <- t(sapply(range.bl, function(targetReturn) {
    w <- MV_QP(x, targetReturn, long_only = long_only)
    c(sd(x %*% w), w)
}))

weigthsl = round(risk.bl[, -1], 4)
colnames(weigthsl) = colnames(x)
weigthsl %>% head()

    DBC IWM    IYR QQQ SPY

[1,] 0.9784 0 0.0016 0 0 [2,] 0.9234 0 0.0566 0 0 [3,] 0.8685 0 0.1115 0 0 [4,] 0.8135 0 0.1665 0 0 [5,] 0.7585 0 0.2215 0 0 [6,] 0.7035 0 0.2765 0 0

risk.bl = risk.bl[, 1]
rets.bl = weigthsl %*% µ.b
fan = 12
plot(x = risk.bl * fan^0.5, y = rets.bl * fan, col = 2, pch = 21, xlab = "Annualized Risk ", 
    ylab = "Annualized Return", main = "long only EF with solve.QP")

Example - Multi Layer Optimization

# demo(multi_layer_optimization)