Tactical Asset Allocation using Machine Learning

Overview

This vignette demonstrates how to construct a machine learning tactical asset allocation strategy using public data and open source. The proposed strategy rotates between two ETFs: high risk such as the SPY ETF and low risk such as 7-10 years IEF Treasury ETF. The strategy returns an improved performance over a naive strategy that allocates 60-40 between the SPY and IEF over time. We provide a detailed description of the implementation, data collection, feature space selection, training, and back-testing. While the literature on return predictability is extensive, it is still important to understand such predictability from the investor’s perspective rather than the econometrician’s alone. We hope that this vignette would encourage further research and reproducibility.

Getting the Data

We refer to the quantmod package to download data from Yahoo Finance. In particular, we focus on a number of market indicators: the SPY ETFs which tracks the S&P 500 index, the VIX index known as the fear gauge which indicates the investors view about the market volatility, the GLD gold ETF, the 7-10 years treasury bond ETF, and XLF the financial sector ETF. Note that all of which are tradable except for the VIX. While VIX can be traded using an exchange traded note, our analysis here does not include its volume. Additionally, the VXX became available in early 2009, the inclusion of which would limit the sample period.

The implementation is conducted mainly using R. In particular, we refer to a number of R packages. We rely on the quantmod package to download market data from Yahoo Finance, the lubridate library to manipulate date formats, plyr for general data manipulation, the glmnet for training and implementation of the machine learning algorithm, the PerformanceAnalytics for financial analytics/summary, and both ggplot2 and plotly for visualizations.

library(quantmod)
library(lubridate)
library(plyr)
library(glmnet)
library(PerformanceAnalytics)
library(ggplot2)
library(plotly)
library(parallel)
rm(list = ls())
t1 <- "1990-01-01"
v <- c("SPY","GLD","IEF","XLF")
P.list <- lapply(v, function(sym) get(getSymbols(sym,from = t1)) )
getSymbols("^VIX",from = t1)

[1] "^VIX"

P.list <- c(P.list,list(VIX))

The P.list returns a list of data for each symbol containing 6 columns. For each item, there are different number of observations. The reason of which is that each ETF/index dates to a different time period.

sapply(P.list,dim)

     [,1] [,2] [,3] [,4] [,5]
[1,] 6836 3861 4443 5346 7615
[2,]    6    6    6    6    6

To see the starting date of each, we can run the following command:

lapply(P.list, function(x)  first(date(x)) )

[[1]]
[1] "1993-01-29"

[[2]]
[1] "2004-11-18"

[[3]]
[1] "2002-07-30"

[[4]]
[1] "1998-12-22"

[[5]]
[1] "1990-01-02"

We note that the SPY ETF has the earliest inception date, whereas the GLD has the latest one.

Feature Space

To construct the main feature space, we focus on the adjusted prices (sixth column) along with volume (fifth column). As mentioned above, note that volume is available for each symbol except the VIX since it is not tradable.

P.list5 <- lapply(P.list, function(x) x[,5])
P.list6 <- lapply(P.list, function(x) x[,6])

Since we are dealing with xts and zoo objects, it is straightforward to merge the time series along:

P5 <- na.omit(Reduce(function(...) merge(...),P.list5 ))
P6 <- na.omit(Reduce(function(...) merge(...),P.list6 ))
# adjust names
names(P5) <- names(P6) <- c("SPY","GLD","IEF","XLF","VIX")
names(P5) <- paste(names(P5),"vol",sep = "_")
summary(P5$VIX_vol)

     Index               VIX_vol 
 Min.   :2004-11-18   Min.   :0  
 1st Qu.:2008-09-19   1st Qu.:0  
 Median :2012-07-19   Median :0  
 Mean   :2012-07-20   Mean   :0  
 3rd Qu.:2016-05-20   3rd Qu.:0  
 Max.   :2020-03-23   Max.   :0  

P5$VIX_vol <- NULL

For each adjusted price, we compute the returns. In addition, we compute the moving average of each return/change over the last 25 days. Our feature space mainly constitutes of the daily return of each symbol, the deviation of each from its MA, and the daily volume of each ETF. Note that R6_roll is a technical analysis tool that demonstrates whether the returns today are too high/low with respect to the corresponding MA. We stack the feature space in the R object. The variable of interest we are trying to model is the next day change in the SPY price

R6 <- Return.calculate(P6)
# add rolling difference 
R6_roll <- R6 - rollapply(R6,25,mean)
names(R6_roll) <- paste(names(R6_roll),"_roll",sep="")
R <- na.omit(merge(R6,R6_roll,P5))
SPY_next <- stats::lag(R$SPY,-1)
names(SPY_next) <- "SPY_next"
R <- na.omit(merge(SPY_next,R))

Finally, the data ranges between

range(date(R))

[1] "2004-12-27" "2020-03-20"

Let’s take a look at the contemporaneous correlation of each feature with the SPY return

cor(R)[,"SPY"]

   SPY_next         SPY         GLD         IEF         XLF         VIX    SPY_roll    GLD_roll    IEF_roll 
-0.11907331  1.00000000  0.02411089 -0.42064985  0.83082759 -0.73426497  0.98736461  0.02890102 -0.41768268 
   XLF_roll    VIX_roll     SPY_vol     GLD_vol     IEF_vol     XLF_vol 
 0.81905191 -0.71892022 -0.13650049 -0.08940766 -0.05839356 -0.05589110 

We observe that the VIX is highly negatively correlated with the same day SPY return. The same holds true for IEF, with -41% correlation. We also note that there is a high positive correlation between the SPY and XLF. Since, the SPY_roll is a function of the same day return, it exhibits high correlation. For the volume, we observe a weak negative correlation. Nonetheless, what matters the more is how the feature space correlates with the next day SPY return than the same return. The following command provides us a perspective on such:

cor(R)[,"SPY_next"]

    SPY_next          SPY          GLD          IEF          XLF          VIX     SPY_roll     GLD_roll 
 1.000000000 -0.119073310 -0.045545477  0.053076235 -0.090127108  0.070958696 -0.117080655 -0.041258439 
    IEF_roll     XLF_roll     VIX_roll      SPY_vol      GLD_vol      IEF_vol      XLF_vol 
 0.048699628 -0.089478329  0.073908555  0.001087906  0.008234174 -0.017657529  0.008429678 

Mainly, the market shows a reversal behavior, where today’s return has a negative 8% correlation with next day. This is also evident for the other indicators. Overall, we witness a weaker correlation between the next day return and the feature space. Obviously, predicting the next day return as a much more challenging task.

Response Variable

Rather than focusing on the next day return, we relate to the change in the SPY price as the response variable with two levels: down (-1) and up (+1). Specifically, we define

R$CHANGE_next <- 1
R$CHANGE_next[R$SPY_next < -0.01] <- -1
table(R$CHANGE)

  -1    1 
 455 3380 

# stack into a dataset rather than an xts object
ds <- data.frame(date = date(R),R)
rownames(ds) <- NULL
ds$SPY_next <- NULL # drop the next day return
# define features
features <- names(ds)[!names(ds) %in% c("date","CHANGE_next")]

Mainly the SPY exhibits 12% of the time a daily drop that is less than -1%. To see whether there is any heterogeneity across the feature space for each level, let’s take a look at the average of the feature space with respect to each level:

sum_change <- dlply(ds,"CHANGE_next",function(x)  apply(x[,features],2,mean) )
(sum_change[[2]] - sum_change[[1]])/abs(sum_change[[1]])

        SPY         GLD         IEF         XLF         VIX    SPY_roll    GLD_roll    IEF_roll    XLF_roll 
 2.11145747 -0.86021267 -0.45423073  0.02854905 -0.18182749 -1.22943200 -1.15739654  1.23701150 -1.19965445 
   VIX_roll     SPY_vol     GLD_vol     IEF_vol     XLF_vol 
 1.19411229 -0.37144701 -0.24487036 -0.09114107 -0.41105495 

Relatively, we observe that the above feature space provides some discrepancy between the two levels. For instance, the previous day return tends to be larger when the market goes up. The opposite is true for gold and bonds.

Machine Learning Application

After defining the feature space and the response variable, we need to predict the probability of the market going up or down on a daily basis. Given this probability we, eventually, will construct our tactical asset allocation strategy.

To get started, we fit a binomial model with an elastic penalty on weekly basis to find the optimal weights to map the feature space into the next day change. Since we use an elastic net, it combines between two penalties. The first is the LASSO which acts as a constraint on the first norm of the weights and serves as an elimination process. The second is the ridge regression which is a constraint on the second norm and serves as a shrinkage approach toward zero. By design, \(\alpha=0.5\) according to the glmnet package acts as in-between penalty of the two approaches. The Lagrangian or the magnitude attributed to the penalty of each constraint, denoted by \(\lambda\), is determined using 10 folds cross validation (henceforth CV).

In the code below, we run a loop in which we train the model using 50 weeks of history and predict the next day change using the following week. Given the 50 weeks of training set, we determine the optimal model using 10 folds CV. Given the test set, we predict the probability of the SPY going either up or down. To avoid data leakage into the test set, we drop the last observation in the training set, as it may contain knowledge about the following week price change.

weeks <- date(unique(floor_date(ds$date,"week")))
weeks <- c(weeks, last(weeks) + weeks(1))
W <- 50
al <- 0.5 # net elastic
ds_predict <- data.frame()
ds_beta <- list()
w_seq <- W:(length(weeks)-2)
ds_predict_f <- function(w) {
  #cat("This is week ",w, " out of ",length(weeks),"\n")
  # training set consists of relatively 250 daily observations
  train.weeks <- weeks[(w-W+1):(w+1)]
  train.index <- which((ds$date > train.weeks[1]) & (ds$date <= train.weeks[W+1]))
  
  # the weekly is around 5 days
  test.weeks <-  weeks[w+1:2]
  test.index <- which((ds$date > test.weeks[1]) & (ds$date < test.weeks[2]))
  
  # drop the last obs from the train set to avoid leakage
  DS <- ds[train.index[-length(train.index)],]
  x_train <- model.matrix( ~ .-1, DS[,features])
  
  # use CV
  set.seed(17)
  try_error <- try(lm <- cv.glmnet(x=x_train,y = as.factor(DS$CHANGE_next), intercept=FALSE,
                 family =   "multinomial", alpha=al, nfolds=10,parallel = T),silent = TRUE)
  
  i <- 1
  while(inherits(try_error,"try-error")) {
    #cat("Error in CV","\n")
    try_error <- try(lm <- cv.glmnet(x=x_train,y = as.factor(DS$CHANGE_next),
                  intercept=FALSE, family =   "multinomial", alpha=al, nfolds=10,
                  parallel = T),silent = TRUE)
    i <- i + 1
    if (i == 10)
        lm <- lm
    }
  
  # assign the lambda
  best_lambda <- lm$lambda.min
  
  # find the optimal model
  lm.star = glmnet(x=x_train,y = as.factor(DS$CHANGE_next), intercept=FALSE ,
                   family =   "multinomial", alpha=al, lambda = best_lambda)
  
  # fit the test sample
  DS_test <- ds[test.index,]
  x_test <- model.matrix( ~ .-1, DS_test[,features])
  DS_predict <- predict(lm.star,x_test,type = "response")
  # stack in data
  DS_predict <- data.frame(DS_predict)
  names(DS_predict) <- c("dn","up")
  DS_predict$date <- ds[test.index,"date"]
  #ds_predict <- rbind(ds_predict,DS_predict)
  
  # finally keep track of the glmnet results in a list
  ds_beta <- c(ds_beta,list(lm.star))
  
  list(DS_predict = DS_predict,lm_list = list(lm.star))
  }
mclapply_list <- mclapply(w_seq,ds_predict_f,mc.cores = detectCores())
ds_predict_l <- lapply(mclapply_list, function(x) x$DS_predict )
ds_predict <- ldply(ds_predict_l,data.frame)
head(ds_predict)

The above code stacks the predicted probabilities in a data frame object named ds_predict. From each iteration, we keep track of the optimal model weights, which allows us to investigate the important features over time. Due to the rolling window nature, which is not recursive in this case, we can refer to the mclapply function for parallelization

summary(ds_predict)

       dn                  up                 date           
 Min.   :0.0000006   Min.   :0.0000156   Min.   :2005-12-12  
 1st Qu.:0.0506136   1st Qu.:0.7717289   1st Qu.:2009-07-08  
 Median :0.1213678   Median :0.8786322   Median :2013-01-31  
 Mean   :0.1558049   Mean   :0.8441951   Mean   :2013-01-30  
 3rd Qu.:0.2282711   3rd Qu.:0.9493864   3rd Qu.:2016-08-24  
 Max.   :0.9999844   Max.   :0.9999994   Max.   :2020-03-20  

Signal Extraction

We merge the predicted probabilities with the original data to construct the main datasheet, which we will refer to for the trading strategy and, hence, the back-testing.

ds2 <- merge(ds,ds_predict,by = "date")
updn <- data.frame(dn = ds2$dn,up =ds2$up)
rownames(updn) <- ds2$date
updn <- as.xts(updn)
updn_roll <- na.omit(rollapply(updn,25,mean))
names(updn_roll) <- c("dn_roll","up_roll")
ds.updn <- data.frame(date = date(updn_roll),updn_roll)
rownames(ds.updn) <- NULL
ds3 <- merge(ds2,ds.updn, by = "date")

To create a smooth signal, we consider the 25 days moving average of the predicted probability. We merge the smoothed probabilities along with the ETF returns in a new data frame named ds3.

To again an initial perspective on the extracted signal, we plot the cumulative return on the SPY along with the probability going down over time. Interestingly, we observe that the signal (red line), which is data driven, reaches the highest level during the 2007-09 financial crisis, during which the SPY (black line) exhibits the largest over the whole sample. Additionally, while the probability seems declining the last couple of year, we observe that it picks up during the recent market turmoils in late 2018. Overall, we observe a negative relationship between the two.

plot(cumsum(SPY) ~ date, data = ds3, type = "l",ylab = "", main  = "cumulative return of SPY versus probability down")
lines(dn_roll ~ date, data = ds3,col = 2)

abline(h = 0, lty=2)
legend("topleft",c("SPY","Probability Down"),col = 1:2, lty = 1)

grid(10)

The above suggests that the probability of down/up can be utilize to avoid significant market drops. The question is how can we utilize this? In the following section, we propose an automated trading strategy that tries to materialize the above insights.

Tactical Asset Allocation Strategy

In this section, we discuss how we implement the tactical allocation strategy using the above extracted signal. In order to demonstrate this, we refer to some mathematical notation. Let \(\hat{\pi}_{t+1}\) denote the forecasted probability that the SPY at \(t+1\) will go up. The hat (\(\hat{\cdot}\)) implies that the probability is a forecast established at time \(t\) rather than \(t+1\). Depending on the intensity of signal, we create a portfolio that either longs the SPY or the IEF. To put formally, the return of the strategy at time \(t+1\) is given by \[\begin{equation} r_{s,t+1} = I_{\left[ \hat{\pi}_{t+1} > a \right]} r^{e}_{t+1} + I_{\left[ \hat{\pi}_{t+1} \leq a \right]} r^{b}_{t+1} \end{equation}\]

where \(r^{e}_{t+1}\) and \(r^{b}_{t+1}\) denote the return on SPY and IEF at \(t+1\), respectively. The parameter \(a\) is a predetermined minimum level of confidence for the investor to participate in the equity market. The strategy implies that the investor goes long in the equity alone, if and only if the confidence level that the market will go up the next day is high enough. If not, the strategy goes $1 in the IEF ETF alone.

Backtesting

As a benchmark, we compare our results with respect to a 60-40 strategy that invests 60% in SPY and 40% in IEF. While the benchmark is re-balanced on a daily basis, it does not incorporate any market information. In other words, regardless of what happens in the market, the benchmark will always maintain a 60-40 allocation. Additionally, we compare the strategy with each ETF alone. As a performance summary, we create a summary function that takes \(a\) as the main input and reports the performance of the strategy with respect to each.

next_f <- function(x) c(x[-1],NA)
ds3$IEF_next <- next_f(ds3$IEF)
ds3$SPY_next <- next_f(ds3$SPY)
ds3$BENCHMARK <- with(ds3,0.6*SPY_next + 0.4*IEF_next)
plot_performance <- function(a) {
  
  # round the probabilities to nearest 0.05
  ds3$up_roll <- round(ds3$up_roll/0.05)*0.05
  ds3$dn_roll <- round(ds3$dn_roll/0.05)*0.05
  ds3$PORT <- with(ds3, (up_roll >= a)*(SPY_next)  +  (up_roll < a)*IEF_next )
  
  
  # load ds3 into a ggplot friendly data
  ds_plot <- data.frame(Date = ds3$date, return = cumsum(ds3$PORT), Type = "Strategy")
  ds_plot <- rbind(ds_plot, 
                   data.frame(Date = ds3$date, return = cumsum(ds3$BENCHMARK), Type = "Benchmark"))
  ds_plot <- rbind(ds_plot, 
                   data.frame(Date = ds3$date, return = cumsum(ds3$SPY_next), Type = "SPY"))
  ds_plot <- rbind(ds_plot, 
                   data.frame(Date = ds3$date, return = cumsum(ds3$IEF_next), Type = "IEF"))
  ds_plot <- rbind(ds_plot, 
                   data.frame(Date = ds3$date, return = ds3$dn_roll, Type = "Probability Down"))
  
  ds_plot <- na.omit(ds_plot)
  
  p <- ggplot(ds_plot) + geom_line(aes(x = Date,y = return,colour = Type)) 
  p <- p + geom_abline(intercept =0,linetype = "dashed")
  p <- ggplotly(p,height =  500, width = 900)
  
  # also return the data for performance comparison
  ds_perf <- ds3[,c("IEF_next","SPY_next","BENCHMARK","PORT","SPY_next")]
  rownames(ds_perf) <- ds3$date
  ds_perf <- as.xts(ds_perf)
  list(plot_perf = p,data_perf = ds_perf)
}

Illustration

Let’s take a look the performance over time when \(a=90\%\). This represents a situation in which the investor allocates his wealth mainly to equity if the confidence level that the SPY is going up the following day is at least 90%. Otherwise, the investor is conservative and allocates his wealth to the Treasury ETF mainly.

perf1 <- plot_performance(0.9)
perf1$plot_perf

In the above plot, we observe that the proposed strategy dominates the other candidates in terms of cumulative return over the period, while at the same time trailing the SPY. Moreover, we observe that the strategy avoids the market crash during the 2007-09 financial crisis, indicating that it has a lower downside risk than the benchmark as well as the SPY. Additionally, given the market turmoil in late 2018, we observe that the the proposed strategy is successful in executing a flight to quality during times of uncertainty.

Suppose that the investor is more conservative, such he only invests in the SPY if the probability is at least 95%. In this case, we have

perf2 <- plot_performance(0.95)
perf2$plot_perf

We note that the proposed strategy avoids both the 2007-09 market crash and the recent market drop in late 2018. However, we observe that the strategy yields almost the same cumulative return as the benchmark.

Suppose that the investor is more risk tolerant, such he is willing to invest in the SPY with a level of confidence of 85%. In this case, we observe that the strategy does well in terms of cumulative return.

perf3 <- plot_performance(0.85)
perf3$plot_perf

However, at the same time, we observe that the strategy fails to shift to Treasury bonds during the recent market turmoils, draining the performance of the strategy over time.

Risk-Adjusted Returns

In terms of risk adjusted returns, we consider a number of statistics. The first is the Jensen’s alpha that captures the abnormal return above the benchmark. The second is Sortino ratio, which adjusts the mean return of the strategy with respect to a downside risk (standard deviation). Additionally, we report the beta of each strategy.

P1 <- perf1$data_perf$PORT
P2 <- perf2$data_perf$PORT
P3 <- perf3$data_perf$PORT
B <- perf1$data_perf$BENCHMARK
B2 <- perf1$data_perf$SPY_next
port_all <- list(P2,P1,P3,B,B2)
port_all <- na.omit(Reduce(merge,port_all))
names(port_all)[1:3] <- c("Portfolio_95","Portfolio_90","Portfolio_85")
table.CAPM(port_all[,1:3],port_all[,4])[c(6,2),]

In the above table, we observe that the more conservative strategies yield an annual alpha of 8%, at the same time, while being beta-neutral On the other hand, we note that when \(a\) decreases, the alpha decreases, and the strategy becomes less market beta-neutral.

In terms of downside risk, we compute the Sortino’s ratio for each strategy. Similar to the Sharpe ratio, we estimate the Sortino’s ratio as the ratio between the mean return and the downside risk. By downside risk, we refer to the semi-standard deviation in which we eliminate the positive returns in calculating the volatility. Finally, we scale the ratio into an annual basis using a scale of \(\sqrt{252}\). As an additional perspective, we report the Sharpe-ratio for comparison:

Sort_ratio <- function(x) sqrt(252)*mean(x)/sd(x[x < 0])
SR <-  function(x) sqrt(252)*mean(x)/sd(x)
Sortino <- apply(port_all,2,Sort_ratio)
Sharpe <- apply(port_all,2,SR)
data.frame(rbind(Sortino,Sharpe))

In all cases, except for the 85% level confidence, we observe that the proposed strategy yields a higher risk-adjusted return than the benchmark. This is more evident for the Sortino than the Sharpe. Hence, this indicates that the strategy outperforms the benchmark by bearing less downside risk, while, at the same time, achieving a better or comparable return.

Summary

This vignette provides a simple cost-efficient trading strategy that deploys machine learning using public data and open source software. While the performance of the strategy depends on a couple of specifications, such as the level of confidence \(a\) or what determines a price change in the market, the strategy is mainly data-driven. One may consider other approaches using cross-validation to determine these inputs automatically. Nonetheless, a similar strategy can be also deployed to screen different stocks or ETFs. For instance, one may consider a similar approach to perform as tactical asset allocation across sector ETFs. i.e. sector rotation strategy. We leave this for future investigation.

Disclaimer

The above vigentte is a short illustration of an ongoing research co-authored with Kris Boudt (Vrije Universiteit Brussel, University of Amsterdam and Finvex), Muzafer Cela (Vrije Universiteit Brussel), and Majeed Simaan (Stevens Institute of Technology). The research is titled “In Search of Return Predictability: Evidence from Machine Learning and Tactical Allocation”. For further application, see the implication of the this analysis into constructing machine learning frontiers via this paper.