A Vector Autoregression Trading Model

The vector autoregression (VAR) framework is common in econometrics for modelling correlated variables with bi-directional relationships and feedback loops. If you google “vector autoregression” you’ll find all sorts of academic papers related to modelling the effects of monetary and fiscal policy on various aspects of the economy. This is only of passing interest to traders.

However, if we consider that the VAR framework finds application in the modelling of correlated time series, the implication being that correlation implies a level of forecasting utility, then perhaps we could model a group of related financial instruments and make predictions that we can translate into trading decisions?

So we’ll give that a try. But first, a brief overview of VAR models.

Overview of VAR models

The univariate autoregression (AR) is a model of a time series as a function of past values of itself:

(Y_t = \alpha + \beta_1 Y_{t-1}+ \beta_2 Y_{t-2} )

That’s an AR(2) model because it uses two previous values in the time series (Y) to estimate the next value. The name of the game is figuring out how many previous values to use, and estimating the coefficients (the (\beta)s) and the intercept ((\alpha)).

A vector autoregression (VAR) is an extension of this idea. It models multiple time series that affect one another together, as a system. It specifically allows for bi-directional relationships such as feedback loops, where say an increase in variable (X) may predict an increase in variable (Y), but equally an increase in variable (Y) may predict an increase in variable (X).

Here’s a VAR(1) model of two time series, (Y_1) and (Y_2):

(Y_{1,t} = \alpha_1 + \beta_{11,1} Y_{1,t-1} + \beta_{12, 1} Y_{2, t-1} + \epsilon_{1, t})
(Y_{2,t} = \alpha_2 + \beta_{21,1} Y_{1,t-1} + \beta_{22, 1} Y_{2, t-1} + \epsilon_{2, t})

The model uses a single lag of each time series to predict the next values of both time series. It requires the estimation of four coefficients and two intercepts.

Just looking at the single lag case, you get a sense that these models have lots of parameters. Which of course triggers all the usual alarm bells around overfitting. In fact, if we have (N) time series and (p) lags in a VAR model, we must estimate (N + pN^2) parameters!

How do you figure out the number of lags?

Standard practice in econometrics is to use an information criterion. It’s questionable how useful that would be in modelling financial asset returns, and in my view it makes sense to stick with a single lag unless you have a compelling reason to do otherwise.

If you must, lean towards the Bayesian information criterion (BIC), which introduces a penalty term for the number of parameters in the model (the Aikake information criterion does too, but the BIC’s penalty is bigger).

A trading problem

Say we have a basket of stocks that we believe to be related in some way. Also, suppose that their relationship implies a degree of predictability among basket constituents. If we could forecast the relative returns between basket constituents, we might have the makings of a trading strategy.

In this post, we’re going to focus on solving the forecasting problem using the VAR framework. In reality, the problem of universe selection (identifying a suitable basket) is a much bigger problem, but we’re going to assume this problem is solved for the sake of the exercise, and that we have a basket of stocks worthy of considering under our VAR framework.

In fact, this is something we’ve solved – in our Machine Learning and Big Data Bootcamp. Go here to find out more and sign up to the wait list for the next edition of Bootcamp.

In this example, our group of stocks appeared in the network model of stock relationships that we built using the Graphical Lasso. This is only a single input into the universe selection model that we trade with, but it will do fine for demonstrating this VAR model. We’ll take the stocks in the little purple cluster consisting of residential construction stocks:

# Extract basket prices
tickers <- c('KBH', 'LEN', 'PHM', 'DHI', 'TOL', 'MTH', 'MDC')

This group is a fairly arbitrary choice – the basket is small enough that we can explore VAR models efficiently but other than that there’s nothing particularly special about it. Other than the fact that the Graphical Lasso identified relationships among the group’s stocks.

You can get historical prices and volumes for these tickers via tidyquant::tq_get, which wraps quantmod::getSymbols:

### VAR MODEL STRATEGY ###

library(tidyquant)
library(tidyverse)

# Load basket prices
tickers <- c('KBH', 'LEN', 'PHM', 'DHI', 'TOL', 'MTH', 'MDC')
basket_prices <- tq_get(tickers, get='stock.prices', from='2000-01-01', to = '2020-01-01') %>%
  rename(ticker = symbol)

basket_prices %>%
  ggplot(aes(x = date, y = adjusted)) +
    geom_line(aes(color = ticker))

Plotting the price series:

 

We’re going to work with returns for our VAR model, so let’s make those now:

# Calculate returns
ret <- basket_prices %>%
  group_by(ticker) %>%
  filter(date >= '2000-01-01', date < '2020-01-01') %>%
  tq_transmute(select=adjusted, mutate_fun=periodReturn, period="daily", type="log")

# make wide dataframe of returns
ret_wide <- ret %>%
  spread(key=ticker, value=daily.returns)

Now we’re ready to build a VAR model. The vars package will do all the heavy lifting for us. So let’s load that, and then construct a VAR(1) model on the first 250 days in our historical data set:

# VAR model
library(vars)

wdw <- 250
var <- VAR(ret_wide[1:wdw, -1], p=1, type="const")

That’s it! Pretty amazing that we can do so much with so little code.

To get a sense for how much work has just been done under the hood, run the command summary(var) and you’ll get a whole lot of output about the model. You’ll see model estimation statistics for each coefficient for each stock in our basket, including insight on what R thought were the significant variables in each equation, as well as information on the residuals of each model.

So now that we’ve got our model, how can we use it to forecast prices? How might we use that information to make trading decisions?

We can get the single step ahead prediction from the model by calling predict, which returns a list with a whole bunch of stuff in addition to the actual forecasts. To extract the forecasts into their own object, we need to do some typically funky R list manipulation:

# step ahead predictions
p <- predict(var, n.ahead = 1) 
fcsts <- t(do.call(rbind, lapply(p$fcst, "[", 1)))

The fcsts object holds a return forecast for each stock in our basket. It looks like this:

> fcsts
              DHI        KBH         LEN         MDC       MTH        PHM          TOL
[1,] -0.007694804 0.01069015 0.008959145 0.008029125 0.0200517 0.01437472 -0.001169301

That’s cool. But how does it compare to the actual return values for the next period?

# predictions vs actuals
ret_wide[wdw+1, -1] - fcsts

# OUTPUT:
#         DHI        KBH        LEN         MDC         MTH        PHM         TOL
# 1 0.04469033 0.03376181 0.04917655 0.007311673 -0.02175964 0.01846445 -0.00187468

Not amazing by the looks of that!

How did we go predicting the correct sign?

# sign of predictions vs actuals
sign(ret_wide[wdw+1, -1]/fcsts)

# OUTPUT
#  DHI KBH LEN MDC MTH PHM TOL
#  1  -1   1   1   1  -1   1   1

Here, a positive number indicates that the return sign was predicted correctly; a negative number indicates an incorrect prediction.

We’ve done a little better here. The model predicted six of eight return signs correctly.

Let’s look at whether we were able to correctly predict the rank of the next period’s returns:

# rank of predictions vs actuals
rank(t(ret_wide[wdw+1, -1]))
rank(fcsts)

# OUTPUT
# > rank(t(ret_wide[wdw+1, -1]))
# 5 6 7 3 2 4 1
# > rank(fcsts)
# 1 5 4 3 7 6 2

Not terrible. In fact, six of eight are accurate enough that if we traded these ranks long-short we’d have done OK on this prediction.

This is all well and good, but so far we’ve only looked at a single prediction. That tells us next to nothing about how useful this approach might be in a trading strategy. Time for some backtesting.

Backtesting a VAR strategy

There are plenty of ways we could act on the predictions of our VAR model. The most familiar of what we looked at above would be trading long-short on the basis of the predicted ranks of our forecast returns.

Here we’re going to do something slightly different that I read about in one of Ernie Chan’s books: we’ll dollar weight assets based on the normalised difference of each asset’s forecast return, and the mean return of the basket:

# dollar weight assets based on normalised difference to forecast mean basket return
dollar_wgts <- (fcsts - mean(fcsts))/sum(abs(fcsts - mean(fcsts)))  # h.t. Ernie Chan
dollar_wgts

# OUTPUT
# DHI        KBH        LEN         MDC       MTH       PHM        TOL
# -0.3177601 0.06405476 0.02810558 0.008791111 0.2584734 0.1405751 -0.1822399

Then, we can calculate a portfolio return by multiplying the actual returns by the dollar weights and summing:

port_ret <- sum(dollar_wgts * ret_wide[wdw+1, -1])
port_ret

# OUTPUT
# -0.002409894

To backtest our VAR forecasts, we’ll set up a rolling window of 500 days on which to estimate the model’s parameters, making a single day-ahead prediction before rolling the window forward by a day and repeating the process, and storing the dollar weights of each asset in a list:

# rolling VAR estimation and prediction
wdw <- 500
dollar_wgts <- vector('list', nrow(ret_wide)-wdw)
for(i in (wdw):(nrow(ret_wide)-1)) {  # don't need to get a forecast for the final data point for backtesting purposes
  
  var <- ret_wide %>%
    slice((i-wdw):(i-1)) %>%
    dplyr::select(-date) %>%
    VAR(p=1, type="const")
  
  p <- var %>% 
    predict(n.ahead = 1)
  
  fcsts <- t(do.call(rbind, lapply(p$fcst, "[", 1)))
  dollar_wgts[[i]] <- (fcsts - mean(fcsts))/sum(abs(fcsts - mean(fcsts))) 
}

Then it’s a matter of extracting the weights from the list, calculating a time series of the cumulative portfolio return and plotting it:

wgts <- bind_rows(lapply(dollar_wgts, as.data.frame))
port_ret <- data.frame('date' = ret_wide[(wdw+1):(nrow(ret_wide)), 'date'], 
                       'ret' = rowSums(wgts * ret_wide[(wdw+1):(nrow(ret_wide)), -1]))  # first wgt is from idx==wdw, assign return from wdw:(wdw+1)

port_ret <- port_ret %>%
  mutate('cum_ret' = cumsum(ret))

port_ret %>%
  ggplot(aes(x = date, y = cum_ret)) +
    geom_line() +
    labs(
      title = 'Cumulative Returns to VAR Trading Model',
      x = 'Date',
      y = 'Return'
    )

Here’s the output:

Costs aren’t considered on this equity curve, and they’d add up to plenty as we’re rebalancing a basket of seven assets on a daily basis. Also bear in mind that the in-sample period for discovering the basket via the Graphical Lasso was 2010-2019 inclusive.

Finally, I wanted to get an idea of the stability of performance with respect to the length of the estimation window. I simply looped over that backtest for different window lengths and stored the results in a list:

# check stability of window length parameter
wdws <- c(250, 500, 750, 1000)
ports <- vector("list")
for (wdw in wdws) {
  dollar_wgts <- vector('list', nrow(ret_wide)-wdw)
  for(i in (wdw):(nrow(ret_wide)-1)) {  # don't need to get a forecast for the final data point for backtesting purposes
    
    var <- ret_wide %>%
      slice((i-wdw):(i-1)) %>%
      dplyr::select(-date) %>%
      VAR(p=1, type="const")
    
    p <- var %>% 
      predict(n.ahead = 1)
    
    fcsts <- t(do.call(rbind, lapply(p$fcst, "[", 1)))
    dollar_wgts[[i]] <- (fcsts - mean(fcsts))/sum(abs(fcsts - mean(fcsts)))  # h.t. Ernie Chan
    
  }
  
  wgts <- bind_rows(lapply(dollar_wgts, as.data.frame))
  port_ret <- data.frame('date' = ret_wide[(wdw+1):(nrow(ret_wide)), 'date'], 
                         'ret' = rowSums(wgts * ret_wide[(wdw+1):(nrow(ret_wide)), -1]))  # first wgt is from idx==wdw, assign return from wdw:(wdw+1)
  
  name <- paste0(wdw)
  port_ret <- port_ret %>%
    mutate(!!name := cumsum(ret))
  
  ports <- c(ports, list(port_ret))
}

With a bit of tidy R manipulation, we can plot the equity curves for the different window lengths:

# plot returns by window length
ports %>%
  map(~ (.x %>% dplyr::select(-ret))) %>%  # drop ret column from each dataframe in list
  reduce(left_join, by="date") %>%  
  pivot_longer(-date, names_to = "window", values_to = "return") %>%
  ggplot(aes(x = date, y = return)) +
    geom_line(aes(color = window)) +
    labs(
      title = 'Return to VAR models of different data windows',
      x = "Date",
      y = "Return"
    )

Here’s the output:

Conclusion

Thus concludes our whirlwind tour of VAR models and their potential use in trading strategies.

In my experience, the problem of finding a basket of related assets that can be exploited in convergence trades like this one is a bigger problem than coming up with the actual trading strategy itself.

And to be fair, there are simpler ways to go about the actual trading than using a VAR model.

However, the VAR model does demonstrate some utility here, and I think it could be used as one of several inputs into a larger ensemble of predictions, rather than a standalone trading strategy.

What do you think? Is the VAR framework attractive to you as a trading tool? Tell me what you think in the comments.