Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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()

plot of chunk unnamed-chunk-4

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()

plot of chunk unnamed-chunk-6

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()

plot of chunk unnamed-chunk-7

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)

plot of chunk unnamed-chunk-9

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))

plot of chunk unnamed-chunk-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 )

plot of chunk unnamed-chunk-13

### 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

plot of chunk unnamed-chunk-13

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) 

plot of chunk unnamed-chunk-13

  ### 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

plot of chunk unnamed-chunk-14

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")

plot of chunk unnamed-chunk-15

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

plot of chunk unnamed-chunk-15

chart.StackedBar(x$BOP.Weight)

plot of chunk unnamed-chunk-15

chart.StackedBar(x$BOP.Value)

plot of chunk unnamed-chunk-15

# 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"

plot of chunk unnamed-chunk-15

# 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")

plot of chunk unnamed-chunk-15

# `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

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")

plot of chunk unnamed-chunk-17

Example - Multi Layer Optimization

# demo(multi_layer_optimization)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.
You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.