1 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)
setwd("~/Documents/GitHub/rob_gordon_2018")

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

symbol	date	open	high	low	close	volume	adjusted
SPY	2017-11-15	256.62	257.22	255.63	256.44	80811500	251.8421
SPY	2017-11-16	257.52	259.04	257.47	258.62	67777000	253.9830
SPY	2017-11-17	258.22	258.59	257.77	257.86	75756800	253.2366
SPY	2017-11-20	258.14	258.52	257.86	258.30	48075500	253.6688
SPY	2017-11-21	259.18	260.20	258.26	259.99	69176800	255.3284
SPY	2017-11-22	260.00	260.15	259.57	259.76	45033400	255.1026
SPY	2017-11-24	260.32	260.48	260.16	260.36	27856500	255.6918
SPY	2017-11-27	260.41	260.75	260.00	260.23	52274900	255.5642
SPY	2017-11-28	260.76	262.90	260.65	262.87	98971700	258.1568
SPY	2017-11-29	263.02	263.63	262.20	262.71	77512100	257.9997

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

4 Calculate Returns

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

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

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

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

5 Growth of $1000 (6mo)

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

symbol	returns	final_value
QQQ	-0.0272053	1090.9494
SPY	-0.0015889	1072.8945
IWM	-0.0038659	1038.9866
IYR	0.0266291	1016.4514
DBC	-0.0536557	996.2755
GLD	-0.0044290	944.2385
EFA	0.0036824	941.7529
TLT	0.0069282	925.4958
EEM	0.0196629	893.9512

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

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

rank	symbol	final_value
1	QQQ	1090.9494
2	SPY	1072.8945
3	IWM	1038.9866
4	IYR	1016.4514
5	DBC	996.2755

6 Portfolio Optimization

6.1 Create Time-Series Object for Optimization

# convert to ts object
stock_returns <- mo_returns %>% 
  # filter(symbol %in% c("SPY", "GLD", "QQQ", "DBC")) %>%
  spread(key = symbol, value = returns) %>%
  tk_xts(date_var = date) %>% 
  na.omit()

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

rm(portf_minvar)
# Create portfolio object
portf_minvar <- portfolio.spec(assets = colnames(stock_returns))
# portf_minvar <- portfolio.spec(assets = colnames(returns.portfolio))
# Add full investment constraint to the portfolio object
portf_minvar <- add.constraint(portfolio=portf_minvar, type="full_investment")
# Add objective to minimize variance
portf_minvar <- add.objective(portfolio=portf_minvar, type="risk", name="var")
# Add long only constraints
portf_minvar <- add.constraint(portfolio=portf_minvar, type="box", min=0, max=1)
 # Run the optimization
opt_gmv <- optimize.portfolio(R=stock_returns, portfolio=portf_minvar, optimize_method="ROI", trace=TRUE)
print(opt_gmv, digits = 5)

PortfolioAnalytics Optimization ***********************************

Call: optimize.portfolio(R = stock_returns, portfolio = portf_minvar, optimize_method = “ROI”, trace = TRUE)

Optimal Weights: DBC EEM EFA GLD IWM IYR QQQ SPY TLT 0.00000 0.00000 0.00000 0.48752 0.00000 0.00000 0.00000 0.13548 0.37700

Objective Measure: StdDev 0.012616

#summary(opt_gmv)
print("Optimal Weights")

[1] “Optimal Weights”

wts <- opt_gmv['weights'] %>% as.data.frame() %>% 
  mutate(.,symbol = rownames(.)) %>% select(symbol, weights)
wts

symbol	weights
DBC	0.0000000
EEM	0.0000000
EFA	0.0000000
GLD	0.4875174
IWM	0.0000000
IYR	0.0000000
QQQ	0.0000000
SPY	0.1354820
TLT	0.3770005

7 Visualize Portfolio Optimization

8 BACKTEST

8.1 Get Historical Stock Returns (3 Years)

hist_returns <- tq_get(stocks
                     , get  = "stock.prices"
                     , from = Sys.Date() - months(36)
                     , to = Sys.Date()) %>% 
  group_by(symbol) %>%
  tq_transmute(select = adjusted, 
                 mutate_fun = periodReturn, 
                 period     = "monthly", 
                 col_rename = "returns") %>% 
  spread(key = symbol, value = returns) %>%
  tk_xts(date_var = date) %>% 
  na.omit()

8.2 Perform Backtest

# https://rdrr.io/cran/PortfolioAnalytics/f/inst/doc/ROI_vignette.pdf
bt_gmv <- optimize.portfolio.rebalancing(R=hist_returns, portfolio=portf_minvar,
                                         optimize_method="ROI",
                                         rebalance_on="months",
                                         training_period=36-6)

print(bt_gmv)

PortfolioAnalytics Optimization with Rebalancing **************************************************

Call: optimize.portfolio.rebalancing(R = hist_returns, portfolio = portf_minvar, optimize_method = “ROI”, rebalance_on = “months”, training_period = 36 - 6)

Number of rebalancing dates: 8 First rebalance date: [1] “2018-04-30” Last rebalance date: [1] “2018-11-14”

Annualized Portfolio Rebalancing Return: [1] -0.02843635

Annualized Portfolio Standard Deviation: [1] 0.08703119

9 Get Performance

9.1 descriptive statistics

symbol	ArithmeticMean	GeometricMean	Kurtosis	LCLMean(0.95)	Maximum	Median	Minimum	Observations	Quartile1	Quartile3	SEMean	Skewness	Stdev	UCLMean(0.95)	Variance
SPY	0.0060	0.0054	0.0073	-0.0146	0.0564	0.0059	-0.0691	13	-0.0016	0.0319	0.0094	-0.7344	0.0340	0.0265	0.0012
IYR	0.0017	0.0013	-0.4254	-0.0176	0.0406	0.0023	-0.0666	13	-0.0239	0.0266	0.0089	-0.6159	0.0319	0.0210	0.0010
QQQ	0.0077	0.0067	-0.0764	-0.0197	0.0876	0.0060	-0.0860	13	-0.0129	0.0280	0.0126	-0.2076	0.0454	0.0351	0.0021
TLT	-0.0057	-0.0059	-1.5196	-0.0190	0.0286	-0.0114	-0.0326	13	-0.0286	0.0131	0.0061	0.1430	0.0220	0.0076	0.0005
GLD	-0.0042	-0.0044	-0.7278	-0.0162	0.0323	-0.0066	-0.0361	13	-0.0208	0.0063	0.0055	0.3733	0.0199	0.0078	0.0004
IWM	0.0039	0.0029	1.3954	-0.0230	0.0616	0.0098	-0.1099	13	-0.0040	0.0256	0.0124	-1.1468	0.0446	0.0308	0.0020
EEM	-0.0076	-0.0086	-0.4608	-0.0354	0.0830	-0.0058	-0.0876	13	-0.0377	0.0197	0.0128	0.1761	0.0461	0.0202	0.0021
DBC	0.0002	-0.0003	-1.1508	-0.0198	0.0342	0.0075	-0.0562	13	-0.0243	0.0278	0.0092	-0.5503	0.0331	0.0202	0.0011
EFA	-0.0040	-0.0046	0.2145	-0.0249	0.0502	0.0037	-0.0813	13	-0.0189	0.0152	0.0096	-0.7049	0.0345	0.0168	0.0012

9.2 annualized returns

symbol	AnnualizedReturn	AnnualizedSharpe(Rf=0%)	AnnualizedStdDev
SPY	0.0671	0.5690	0.1179
IYR	0.0152	0.1371	0.1107
QQQ	0.0837	0.5325	0.1571
TLT	-0.0690	-0.9040	0.0763
GLD	-0.0516	-0.7486	0.0689
IWM	0.0359	0.2328	0.1543
EEM	-0.0983	-0.6160	0.1596
DBC	-0.0034	-0.0300	0.1146
EFA	-0.0539	-0.4514	0.1194

9.3 downside risks

symbol	DownsideDeviation(0%)	DownsideDeviation(MAR=10%)	DownsideDeviation(Rf=0%)	GainDeviation	HistoricalES(95%)	HistoricalVaR(95%)	LossDeviation	MaximumDrawdown	ModifiedES(95%)	ModifiedVaR(95%)	SemiDeviation
SPY	0.0230	0.0269	0.0230	0.0178	-0.0691	-0.0495	0.0279	0.0706	-0.0673	-0.0543	0.0257
IYR	0.0228	0.0273	0.0228	0.0146	-0.0666	-0.0448	0.0235	0.0959	-0.0637	-0.0542	0.0237
QQQ	0.0277	0.0318	0.0277	0.0303	-0.0860	-0.0589	0.0324	0.1133	-0.0861	-0.0666	0.0315
TLT	0.0185	0.0244	0.0185	0.0085	-0.0326	-0.0313	0.0084	0.0868	-0.0430	-0.0403	0.0145
GLD	0.0152	0.0213	0.0152	0.0107	-0.0361	-0.0279	0.0109	0.1166	-0.0382	-0.0339	0.0124
IWM	0.0330	0.0367	0.0330	0.0217	-0.1099	-0.0670	0.0438	0.1339	-0.1043	-0.0782	0.0346
EEM	0.0353	0.0406	0.0353	0.0282	-0.0876	-0.0704	0.0263	0.2276	-0.0919	-0.0786	0.0306
DBC	0.0246	0.0295	0.0246	0.0117	-0.0562	-0.0547	0.0171	0.1098	-0.0613	-0.0576	0.0247
EFA	0.0279	0.0326	0.0279	0.0154	-0.0813	-0.0615	0.0275	0.1373	-0.0803	-0.0647	0.0258

9.4 downside risk ratios

symbol	Annualiseddownsiderisk	Downsidepotential	Monthlydownsiderisk	Omega	Omega-sharperatio	Sortinoratio	Upsidepotential	Upsidepotentialratio
SPY	0.0795	0.0103	0.0230	1.5770	0.5770	0.2600	0.0163	0.5693
IYR	0.0788	0.0116	0.0228	1.1497	0.1497	0.0761	0.0133	0.6734
QQQ	0.0959	0.0131	0.0277	1.5875	0.5875	0.2769	0.0207	0.7542
TLT	0.0640	0.0129	0.0185	0.5567	-0.4433	-0.3093	0.0072	0.6175
GLD	0.0526	0.0104	0.0152	0.5959	-0.4041	-0.2778	0.0062	1.1078
IWM	0.1142	0.0138	0.0330	1.2820	0.2820	0.1181	0.0177	0.5410
EEM	0.1221	0.0223	0.0353	0.6591	-0.3409	-0.2157	0.0147	0.6629
DBC	0.0852	0.0141	0.0246	1.0159	0.0159	0.0091	0.0143	0.5856
EFA	0.0967	0.0150	0.0279	0.7306	-0.2694	-0.1449	0.0110	0.4958

9.5 value at risk

symbol	VaR
SPY	-0.0543204
IYR	-0.0541734
QQQ	-0.0666155
TLT	-0.0403022
GLD	-0.0338613
IWM	-0.0782136
EEM	-0.0785758
DBC	-0.0575836
EFA	-0.0646908

9.6 sharpe ratio

symbol	ESSharpe(Rf=0%,p=95%)	StdDevSharpe(Rf=0%,p=95%)	VaRSharpe(Rf=0%,p=95%)
SPY	0.0887032	0.1753105	0.1098641
IYR	0.0271849	0.0542334	0.0319830
QQQ	0.0890736	0.1690699	0.1151180
TLT	-0.1327730	-0.2593925	-0.1417623
GLD	-0.1103754	-0.2121981	-0.1246630
IWM	0.0373081	0.0873565	0.0497653
EEM	-0.0827627	-0.1650276	-0.0967537
DBC	0.0036584	0.0067779	0.0038938
EFA	-0.0503964	-0.1173963	-0.0625411

10 Appendix

10.1 Resources & Notes

10.2 Example of Efficiency Frontier from Stack Overflow

library(PortfolioAnalytics)
  returns <- stock_returns
#  define moment function
  annualized.moments <- function(R, scale=12, portfolio=NULL){
    out <- list()
    out$mu <-    matrix(colMeans(R), ncol=1)
    out$sigma <- cov(R)/scale
    return(out)
  }
# define portfolio
  prt <- portfolio.spec(assets=colnames(stock_returns))
  prt <- add.constraint(portfolio=prt, type="long_only")
  #  leverage defaults to weight_sum = 1 so is equivalent to full_investment constraint
  prt <- add.constraint(portfolio=prt, type="leverage")
  prt <- add.objective(portfolio=prt, type="risk", name="StdDev")
# calculate and plot efficient frontier
  prt_ef <- create.EfficientFrontier(R=12*stock_returns, portfolio=prt, type="mean-StdDev", 
                                      match.col = "StdDev", momentFUN="annualized.moments", scale=12)
  xlim <- range(prt_ef$frontier[,2])*c(1, 1.5)
  ylim <- range(prt_ef$frontier[,1])*c(.80, 1.05)
  chart.EfficientFrontier(prt_ef, match.col="StdDev", chart.assets = FALSE, 
                          labels.assets = FALSE, xlim=xlim, ylim=ylim )
  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(returns), cex=.8, pos=4)

  chart.EF.Weights(prt_ef, match.col="StdDev")

10.3 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)),.98,-1.02)
  }
  sol  <- solve.QP(Dmat=Sig, dvec, t(Amat), bvec, meq=meq)$solution 
}

steps=50
x=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    EEM EFA    GLD IWM IYR QQQ SPY    TLT

[1,] 0 0.9395 0 0.0000 0 0 0 0 0.0805 [2,] 0 0.7745 0 0.0000 0 0 0 0 0.2455 [3,] 0 0.6096 0 0.0000 0 0 0 0 0.4104 [4,] 0 0.4446 0 0.0000 0 0 0 0 0.5754 [5,] 0 0.3022 0 0.0285 0 0 0 0 0.6893 [6,] 0 0.2305 0 0.1466 0 0 0 0 0.6429

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

10.4 Example - Multi Layer Optimization

# demo(multi_layer_optimization)

Rob Gordon Rebalancing Tool