An Out of Sample Update on DDN’s Volatility Momentum Trading Strategy and Beta Convexity

The first part of this post is a quick update on Tony Cooper’s of Double Digit Numerics’s volatility ETN momentum strategy from the volatility made simple blog (which has stopped updating as of a year and a half ago). The second part will cover Dr. Jonathan Kinlay’s Beta Convexity concept.

So, now that I have the ability to generate a term structure and constant expiry contracts, I decided to revisit some of the strategies on Volatility Made Simple and see if any of them are any good (long story short: all of the publicly detailed ones aren’t so hot besides mine–they either have a massive drawdown in-sample around the time of the crisis, or a massive drawdown out-of-sample).

Why this strategy? Because it seemed different from most of the usual term structure ratio trades (of which mine is an example), so I thought I’d check out how it did since its first publishing date, and because it’s rather easy to understand.

Here’s the strategy:

Take XIV, VXX, ZIV, VXZ, and SHY (this last one as the “risk free” asset), and at the close, invest in whichever has had the highest 83 day momentum (this was the result of optimization done on volatilityMadeSimple).

Here’s the code to do this in R, using the Quandl EOD database. There are two variants tested–observe the close, buy the close (AKA magical thinking), and observe the close, buy tomorrow’s close.



symbols <- c("XIV", "VXX", "ZIV", "VXZ", "SHY")

prices <- list()
for(i in 1:length(symbols)) {
  price <- Quandl(paste0("EOD/", symbols[i]), start_date="1990-12-31", type = "xts")$Adj_Close
  colnames(price) <- symbols[i]
  prices[[i]] <- price
prices <- na.omit(, prices))
returns <- na.omit(Return.calculate(prices))

# find highest asset, assign column names
topAsset <- function(row, assetNames) {
  out <- row==max(row, na.rm = TRUE)
  names(out) <- assetNames
  out <- data.frame(out)

# compute momentum
momentums <- na.omit(xts(apply(prices, 2, ROC, n = 83),

# find highest asset each day, turn it into an xts
highestMom <- apply(momentums, 1, topAsset, assetNames = colnames(momentums))
highestMom <- xts(t(, highestMom)),

# observe today's close, buy tomorrow's close
buyTomorrow <- na.omit(xts(rowSums(returns * lag(highestMom, 2)),

# observe today's close, buy today's close (aka magic thinking)
magicThinking <- na.omit(xts(rowSums(returns * lag(highestMom)),

out <- na.omit(cbind(buyTomorrow, magicThinking))
colnames(out) <- c("buyTomorrow", "magicalThinking")

# results
charts.PerformanceSummary(out['2014-04-11::'], legend.loc = 'top')
rbind(table.AnnualizedReturns(out['2014-04-11::']), maxDrawdown(out['2014-04-11::']))

Pretty simple.

Here are the results.


> rbind(table.AnnualizedReturns(out['2014-04-11::']), maxDrawdown(out['2014-04-11::']))
                          buyTomorrow magicalThinking
Annualized Return          -0.0320000       0.0378000
Annualized Std Dev          0.5853000       0.5854000
Annualized Sharpe (Rf=0%)  -0.0547000       0.0646000
Worst Drawdown              0.8166912       0.7761655

Looks like this strategy didn’t pan out too well. Just a daily reminder that if you’re using fine grid-search to select a particularly good parameter (EG n = 83 days? Maybe 4 21-day trading months, but even that would have been n = 82), you’re asking for a visit from, in the words of Mr. Tony Cooper, a visit from the grim reaper.


Moving onto another topic, whenever Dr. Jonathan Kinlay posts something that I think I can replicate that I’d be very wise to do so, as he is a very skilled and experienced practitioner (and also includes me on his blogroll).

A topic that Dr. Kinlay covered is the idea of beta convexity–namely, that an asset’s beta to a benchmark may be different when the benchmark is up as compared to when it’s down. Essentially, it’s the idea that we want to weed out firms that are what I’d deem as “losers in disguise”–I.E. those that act fine when times are good (which is when we really don’t care about diversification, since everything is going up anyway), but do nothing during bad times.

The beta convexity is calculated quite simply: it’s the beta of an asset to a benchmark when the benchmark has a positive return, minus the beta of an asset to a benchmark when the benchmark has a negative return, then squaring the difference. That is, (beta_bench_positive – beta_bench_negative) ^ 2.

Here’s some R code to demonstrate this, using IBM vs. the S&P 500 since 1995.

ibm <- Quandl("EOD/IBM", start_date="1995-01-01", type = "xts")
ibmRets <- Return.calculate(ibm$Adj_Close)

spy <- Quandl("EOD/SPY", start_date="1995-01-01", type = "xts")
spyRets <- Return.calculate(spy$Adj_Close)

rets <- na.omit(cbind(ibmRets, spyRets))
colnames(rets) <- c("IBM", "SPY")

betaConvexity <- function(Ra, Rb) {
  positiveBench <- Rb[Rb > 0]
  assetPositiveBench <- Ra[index(positiveBench)]
  positiveBeta <- CAPM.beta(Ra = assetPositiveBench, Rb = positiveBench)
  negativeBench <- Rb[Rb < 0]
  assetNegativeBench <- Ra[index(negativeBench)]
  negativeBeta <- CAPM.beta(Ra = assetNegativeBench, Rb = negativeBench)
  out <- (positiveBeta - negativeBeta) ^ 2

betaConvexity(rets$IBM, rets$SPY)

For the result:

> betaConvexity(rets$IBM, rets$SPY)
[1] 0.004136034

Thanks for reading.

NOTE: I am always looking to network, and am currently actively looking for full-time opportunities which may benefit from my skill set. If you have a position which may benefit from my skills, do not hesitate to reach out to me. My LinkedIn profile can be found here.

Constant Expiry VIX Futures (Using Public Data)

This post will be about creating constant expiry (E.G. a rolling 30-day contract) using VIX settlement data from the CBOE and the spot VIX calculation (from Yahoo finance, or wherever else). Although these may be able to be traded under certain circumstances, this is not always the case (where the desired expiry is shorter than the front month’s time to expiry).

The last time I visited this topic, I created a term structure using publicly available data from the CBOE, along with an external expiry calendar.

The logical next step, of course, is to create constant-expiry contracts, which may or may not be tradable (if your contract expiry is less than 30 days, know that the front month has days in which the time to expiry is more than 30 days).

So here’s where we left off: a way to create a continuous term structure using CBOE settlement VIX data.

So from here, before anything, we need to get VIX data. And while the getSymbols command used to be easier to use, because Yahoo broke its API (what else do you expect from an otherwise-irrelevant, washed-up web 1.0 dinosaur?), it’s not possible to get free Yahoo data at this point in time (in the event that Josh Ulrich doesn’t fix this issue in the near future, I’m open to suggestions for other free sources of data which provide data of reputable quality), so we need to get VIX data from elsewhere (particularly, the CBOE itself, which is a one-stop shop for all VIX-related data…and most likely some other interesting futures as well.)

So here’s how to get VIX data from the CBOE (thanks, all you awesome CBOE people! And a shoutout to all my readers from the CBOE, I’m sure some of you are from there).

VIX <- fread("", skip = 1)
VIXdates <- VIX$Date
VIX$Date <- NULL; VIX <- xts(VIX,, format = '%m/%d/%Y'))
spotVix <- Cl(VIX)

Next, there’s a need for some utility functions to help out with identifying which futures contracts to use for constructing synthetics.

# find column with greatest days to expiry less than or equal to desired days to expiry
shortDurMinMax <- function(row, daysToExpiry) {
  return(max(which(row <= daysToExpiry)))

# find column with least days to expiry greater desired days to expiry
longDurMinMax <- function(row, daysToExpiry) {
  return(min(which(row > daysToExpiry)))

# gets the difference between the two latest contracts (either expiry days or price)
getLastDiff <- function(row) {
  indices <- rev(which(!
  out <- row[indices[1]] - row[indices[2]]

# gets the rightmost non-NA value of a row
getLastValue <- function(row) {
  indices <- rev(which(!
  out <- row[indices[1]]

The first two functions are to determine short-duration and long-duration contracts. Simply, provided a row of data and the desired constant time to expiry, the first function finds the contract with a time closest to expiry less than or equal to the desired amount, while the second function does the inverse.

The next two functions are utilized in the scenario of a function whose time to expiry is greater than the expiry of the longest trading contract. Such a synthetic would obviously not be able to be traded, but can be created for the purposes of using as an indicator. The third function gets the last two non-NA values in a row (I.E. the two last prices, the two last times to expiry), and the fourth one simply gets the rightmost non-NA value in a row.

The algorithm to create a synthetic constant-expiry contract/indicator is divided into three scenarios:

One, in which the desired time to expiry of the contract is shorter than the front month, such as a constant 30-day expiry contract, when the front month has more than 30 days to maturity (such as on Nov 17, 2016), at which point, the weight will be the desired time to expiry over the remaining time to expiry in the front month, and the remainder in spot VIX (another asset that cannot be traded, at least conventionally).

The second scenario is one in which the desired time to expiry is longer than the last traded contract. For instance, if the desire was to create a contract
with a year to expiry when the furthest out is eight months, there obviously won’t be data for such a contract. In such a case, the algorithm is to compute the linear slope between the last two available contracts, and add the extrapolated product of the slope multiplied by the time remaining between the desired and the last contract to the price of the last contract.

Lastly, the third scenario (and the most common one under most use cases) is that of the synthetic for which there is both a trading contract that has less time to expiry than the desired constant rate, and one with more time to expiry. In this instance, a matter of linear algebra (included in the comments) denotes the weight of the short expiry contract, which is (desired – expiry_long)/(expiry_short – expiry_long).

The algorithm iterates through all three scenarios, and due to the mechanics of xts automatically sorting by timestamp, one obtains an xts object in order of dates of a synthetic, constant expiry futures contract.

Here is the code for the function.

constantExpiry <- function(spotVix, termStructure, expiryStructure, daysToExpiry) {
  # Compute synthetics that are too long (more time to expiry than furthest contract)
  # can be Inf if no column contains values greater than daysToExpiry (I.E. expiry is 3000 days)
  suppressWarnings(longCol <- xts(apply(expiryStructure, 1, longDurMinMax, daysToExpiry),
  longCol[longCol == Inf] <- 10
  # xts for too long to expiry -- need a NULL for rbinding if empty
  tooLong <- NULL
  # Extend the last term structure slope an arbitrarily long amount of time for those with too long expiry
  tooLongIdx <- index(longCol[longCol==10])
  if(length(tooLongIdx) > 0) {
    tooLongTermStructure <- termStructure[tooLongIdx]
    tooLongExpiryStructure <- expiryStructure[tooLongIdx]
    # difference in price/expiry for longest two contracts, use it to compute a slope
    priceDiff <- xts(apply(tooLongTermStructure, 1, getLastDiff), = tooLongIdx)
    expiryDiff <- xts(apply(tooLongExpiryStructure, 1, getLastDiff), = tooLongIdx)
    slope <- priceDiff/expiryDiff
    # get longest contract price and compute additional days to expiry from its time to expiry 
    # I.E. if daysToExpiry is 180 and longest is 120, additionalDaysToExpiry is 60
    maxDaysToExpiry <- xts(apply(tooLongExpiryStructure, 1, max, na.rm = TRUE), = tooLongIdx)
    longestContractPrice <- xts(apply(tooLongTermStructure, 1, getLastValue), = tooLongIdx)
    additionalDaysToExpiry <- daysToExpiry - maxDaysToExpiry
    # add slope multiplied by additional days to expiry to longest contract price
    tooLong <- longestContractPrice + additionalDaysToExpiry * slope
  # compute synthetics that are too short (less time to expiry than shortest contract)
  # can be -Inf if no column contains values less than daysToExpiry (I.E. expiry is 5 days)
  suppressWarnings(shortCol <- xts(apply(expiryStructure, 1, shortDurMinMax, daysToExpiry),
  shortCol[shortCol == -Inf] <- 0
  # xts for too short to expiry -- need a NULL for rbinding if empty
  tooShort <- NULL
  tooShortIdx <- index(shortCol[shortCol==0])
  if(length(tooShortIdx) > 0) {
    tooShort <- termStructure[,1] * daysToExpiry/expiryStructure[,1] + spotVix * (1 - daysToExpiry/expiryStructure[,1])
    tooShort <- tooShort[tooShortIdx]
  # compute everything in between (when time to expiry is between longest and shortest)
  # get unique permutations for contracts that term structure can create
  colPermutes <- cbind(shortCol, longCol)
  colnames(colPermutes) <- c("short", "long")
  colPermutes <- colPermutes[colPermutes$short > 0,]
  colPermutes <- colPermutes[colPermutes$long < 10,]
  regularSynthetics <- NULL
  # if we can construct synthetics from regular futures -- someone might enter an extremely long expiry
  # so this may not always be the case
  if(nrow(colPermutes) > 0) {
    # pasting long and short expiries into a single string for easier subsetting
    shortLongPaste <- paste(colPermutes$short, colPermutes$long, sep="_")
    uniqueShortLongPaste <- unique(shortLongPaste)
    regularSynthetics <- list()
    for(i in 1:length(uniqueShortLongPaste)) {
      # get unique permutation of short-expiry and long-expiry contracts
      permuteSlice <- colPermutes[which(shortLongPaste==uniqueShortLongPaste[i]),]
      expirySlice <- expiryStructure[index(permuteSlice)]
      termStructureSlice <- termStructure[index(permuteSlice)]
      # what are the parameters?
      shortCol <- unique(permuteSlice$short); longCol <- unique(permuteSlice$long)
      # computations -- some linear algebra
      # S/L are weights, ex_S/ex_L are time to expiry
      # D is desired constant time to expiry
      # S + L = 1
      # L = 1 - S
      # S + (1-S) = 1
      # ex_S * S + ex_L * (1-S) = D
      # ex_S * S + ex_L - ex_L * S = D
      # ex_S * S - ex_L * S = D - ex_L
      # S(ex_S - ex_L) = D - ex_L
      # S = (D - ex_L)/(ex_S - ex_L)
      weightShort <- (daysToExpiry - expirySlice[, longCol])/(expirySlice[, shortCol] - expirySlice[, longCol])
      weightLong <- 1 - weightShort
      syntheticValue <- termStructureSlice[, shortCol] * weightShort + termStructureSlice[, longCol] * weightLong
      regularSynthetics[[i]] <- syntheticValue
    regularSynthetics <-, regularSynthetics)
  out <- rbind(tooShort, regularSynthetics, tooLong)
  colnames(out) <- paste0("Constant_", daysToExpiry)

And here’s how to use it:

constant30 <- constantExpiry(spotVix = vixSpot, termStructure = termStructure, expiryStructure = expiryStructure, daysToExpiry = 30)
constant180 <- constantExpiry(spotVix = vixSpot, termStructure = termStructure, expiryStructure = expiryStructure, daysToExpiry = 180)

constantTermStructure <- cbind(constant30, constant180)

chart.TimeSeries(constantTermStructure, legend.loc = 'topright', main = "Constant Term Structure")

With the result:

Which means that between the CBOE data itself, and this function that creates constant expiry futures from CBOE spot and futures prices, one can obtain any futures contract, whether real or synthetic, to use as an indicator for volatility trading strategies. This allows for exploration of a wide variety of volatility trading strategies.

Thanks for reading.

NOTE: I am always interested in networking and hearing about full-time opportunities related to my skill set. My linkedin can be found here.

Furthermore, if you are a volatility ETF/futures trading professional, I am interested in possible project-based collaboration. If you are interested, please contact me.

Creating a VIX Futures Term Structure In R From Official CBOE Settlement Data

This post will be detailing a process to create a VIX term structure from freely available CBOE VIX settlement data and a calendar of freely obtainable VIX expiry dates. This has applications for volatility trading strategies.

So this post, as has been the usual for quite some time, will not be about a strategy, but rather, a tool that can be used for exploring future strategies. Particularly, volatility strategies–which seems to have been a hot topic on this blog some time ago (and might very well be still, especially since the Volatility Made Simple blog has just stopped tracking open-sourced strategies for the past year).

This post’s topic is the VIX term structure–that is, creating a set of continuous contracts–properly rolled according to VIX contract specifications, rather than a hodgepodge of generic algorithms as found on some other websites. The idea is, as of the settlement of a previous day (or whenever the CBOE actually releases their data), you can construct a curve of contracts, and see if it’s in contango (front month cheaper than next month and so on) or backwardation (front month more expensive than next month, etc.).

The first (and most code-intensive) part of the procedure is fairly simple–map the contracts to an expiration date, then put their settlement dates and times to expiry into two separate xts objects, with one column for each contract.

The expiries text file is simply a collection of copied and pasted expiry dates from this site. It includes the January 2018 expiration date. Here is what it looks like:

> head(expiries)
  V1       V2   V3
1 18  January 2006
2 15 February 2006
3 22    March 2006
4 19    April 2006
5 17      May 2006
6 21     June 2006

# 06 through 17
years <- c(paste0("0", c(6:9)), as.character(c(10:17)))

# futures months
futMonths <- c("F", "G", "H", "J", "K", "M",
            "N", "Q", "U", "V", "X", "Z")

# expiries come from
expiries <- read.table("expiries.txt", header = FALSE, sep = " ")

# convert expiries into dates in R
dateString <- paste(expiries$V3, expiries$V2, expiries$V1, sep = "-")
dates <- as.Date(dateString, format = "%Y-%B-%d")

# map futures months to numbers for dates
monthMaps <- cbind(futMonths, c("01", "02", "03", "04", "05", "06",
                                   "07", "08", "09", "10", "11", "12"))
monthMaps <- data.frame(monthMaps)
colnames(monthMaps) <- c("futureStem", "monthNum")

dates <- data.frame(dates)
dates$dateMon <- substr(dates$dates, 1, 7)

contracts <- expand.grid(futMonths, years)
contracts <- paste0(contracts[,1], contracts[,2])
contracts <- c(contracts, "F18")
stem <- ""
#contracts <- paste0(stem, contracts, "_VX.csv")

masterlist <- list()
timesToExpiry <- list()
for(i in 1:length(contracts)) {
  # obtain data
  contract <- contracts[i]
  dataFile <- paste0(stem, contract, "_VX.csv")
  expiryYear <- paste0("20",substr(contract, 2, 3))
  expiryMonth <- monthMaps$monthNum[monthMaps$futureStem == substr(contract,1,1)]
  expiryDate <- dates$dates[dates$dateMon == paste(expiryYear, expiryMonth, sep="-")]
  data <- suppressWarnings(fread(dataFile))
  # create dates
  dataDates <- as.Date(data$`Trade Date`, format = '%m/%d/%Y')
  # create time to expiration xts
  toExpiry <- xts(expiryDate - dataDates,
  colnames(toExpiry) <- contract
  timesToExpiry[[i]] <- toExpiry
  # get settlements
  settlement <- xts(data$Settle,
  colnames(settlement) <- contract
  masterlist[[i]] <- settlement

# cbind outputs
masterlist <-, masterlist)
timesToExpiry <-, timesToExpiry)

# NA out zeroes in settlements
masterlist[masterlist==0] <- NA

From there, we need to visualize how many contracts are being traded at once on any given day (I.E. what’s a good steady state number for the term structure)?

sumNonNA <- function(row) {

simultaneousContracts <- xts(apply(masterlist, 1, sumNonNA),

The result looks like this:

So, 8 contracts (give or take) at any given point in time. This is confirmed by the end of the master list of settlements.

> dim(masterlist)
[1] 3002  145
> tail(masterlist[,135:145])
           H17    J17    K17    M17    N17    Q17    U17    V17    X17    Z17   F18
2017-04-18  NA 14.725 14.325 14.525 15.175 15.475 16.225 16.575 16.875 16.925    NA
2017-04-19  NA 14.370 14.575 14.525 15.125 15.425 16.175 16.575 16.875 16.925    NA
2017-04-20  NA     NA 14.325 14.325 14.975 15.375 16.175 16.575 16.875 16.900    NA
2017-04-21  NA     NA 14.325 14.225 14.825 15.175 15.925 16.350 16.725 16.750    NA
2017-04-24  NA     NA 12.675 13.325 14.175 14.725 15.575 16.025 16.375 16.475 17.00
2017-04-25  NA     NA 12.475 13.125 13.975 14.425 15.225 15.675 16.025 16.150 16.75

Using this information, an algorithm can create eight continuous contracts, ranging from front month to eight months out. The algorithm starts at the first day of the master list to the first expiry, then moves between expiry windows, and just appends the front month contract, and the next seven contracts to a list, before rbinding them together, and does the same with the expiry structure.

termStructure <- list()
expiryStructure <- list()
masterDates <- unique(c(first(index(masterlist)), dates$dates[dates$dates %in% index(masterlist)], Sys.Date()-1))
for(i in 1:(length(masterDates)-1)) {
  subsetDates <- masterDates[c(i, i+1)]
  dateRange <- paste(subsetDates[1], subsetDates[2], sep="::")
  subset <- masterlist[dateRange,c(i:(i+7))]
  subset <- subset[-1,]
  expirySubset <- timesToExpiry[index(subset), c(i:(i+7))]
  colnames(subset) <- colnames(expirySubset) <- paste0("C", c(1:8))
  termStructure[[i]] <- subset
  expiryStructure[[i]] <- expirySubset

termStructure <-, termStructure)
expiryStructure <-, expiryStructure)

Again, one more visualization of when we have a suitable number of contracts:

simultaneousContracts <- xts(apply(termStructure, 1, sumNonNA),

And in order to preserve the most data, we’ll cut the burn-in period off when we first have 7 contracts trading at once.

first(index(simultaneousContracts)[simultaneousContracts >= 7])
termStructure <- termStructure["2006-10-23::"]
expiryStructure <- expiryStructure[index(termStructure)]

So there you have it–your continuous VIX futures contract term structure, as given by the official CBOE settlements. While some may try and simulate a trading strategy based on these contracts, I myself prefer to use them as indicators or features to a model that would rather buy XIV or VXX.

One last trick, for those that want to visualize things, a way to actually visualize the term structure on any given day, in particular, the most recent one in the term structure.

plot(t(coredata(last(termStructure))), type = 'b')

A clear display of contango.

A post on how to compute synthetic constant-expiry contracts (EG constant 30 day expiry contracts) will be forthcoming in the near future.

Thanks for reading.

NOTE: I am currently interested in networking and full-time positions which may benefit from my skills. I may be contacted at my LinkedIn profile found here.

How well can you scale your strategy?

This post will deal with a quick, finger in the air way of seeing how well a strategy scales–namely, how sensitive it is to latency between signal and execution, using a simple volatility trading strategy as an example. The signal will be the VIX/VXV ratio trading VXX and XIV, an idea I got from Volatility Made Simple’s amazing blog, particularly this post. The three signals compared will be the “magical thinking” signal (observe the close, buy the close, named from the ruleOrderProc setting in quantstrat), buy on next-day open, and buy on next-day close.

Let’s get started.


         destfile="longVXX.txt") #requires downloader package
getSymbols('^VIX', from = '1990-01-01')

xiv <- xts(read.zoo("longXIV.txt", format="%Y-%m-%d", sep=",", header=TRUE))
vxx <- xts(read.zoo("longVXX.txt", format="%Y-%m-%d", sep=",", header=TRUE))
vxv <- xts(read.zoo("vxvData.csv", header=TRUE, sep=",", format="%m/%d/%Y", skip=2))
vixVxv <- Cl(VIX)/Cl(vxv)

xiv <- xts(read.zoo("longXIV.txt", format="%Y-%m-%d", sep=",", header=TRUE))
vxx <- xts(read.zoo("longVXX.txt", format="%Y-%m-%d", sep=",", header=TRUE))

vxxCloseRets <- Return.calculate(Cl(vxx))
vxxOpenRets <- Return.calculate(Op(vxx))
xivCloseRets <- Return.calculate(Cl(xiv))
xivOpenRets <- Return.calculate(Op(xiv))

vxxSig <- vixVxv > 1
xivSig <- 1-vxxSig

magicThinking <- vxxCloseRets * lag(vxxSig) + xivCloseRets * lag(xivSig)
nextOpen <- vxxOpenRets * lag(vxxSig, 2) + xivOpenRets * lag(xivSig, 2)
nextClose <- vxxCloseRets * lag(vxxSig, 2) + xivCloseRets * lag(xivSig, 2)
tradeWholeDay <- (nextOpen + nextClose)/2

compare <- na.omit(cbind(magicThinking, nextOpen, nextClose, tradeWholeDay))
colnames(compare) <- c("Magic Thinking", "Next Open", 
                       "Next Close", "Execute Through Next Day")
      maxDrawdown(compare), CalmarRatio(compare))

chart.TimeSeries(log(cumprod(1+compare), base = 10), legend.loc='topleft', ylab='log base 10 of additional equity',
                 main = 'VIX vx. VXV different execution times')

So here’s the run-through. In addition to the magical thinking strategy (observe the close, buy that same close), I tested three other variants–a variant which transacts the next open, a variant which transacts the next close, and the average of those two. Effectively, I feel these three could give a sense of a strategy’s performance under more realistic conditions–that is, how well does the strategy perform if transacted throughout the day, assuming you’re managing a sum of money too large to just plow into the market in the closing minutes (and if you hope to get rich off of trading, you will have a larger sum of money than the amount you can apply magical thinking to). Ideally, I’d use VWAP pricing, but as that’s not available for free anywhere I know of, that means that readers can’t replicate it even if I had such data.

In any case, here are the results.

Equity curves:

Log scale (for Mr. Tony Cooper and others):


                          Magic Thinking Next Open Next Close Execute Through Next Day
Annualized Return               0.814100 0.8922000  0.5932000                 0.821900
Annualized Std Dev              0.622800 0.6533000  0.6226000                 0.558100
Annualized Sharpe (Rf=0%)       1.307100 1.3656000  0.9529000                 1.472600
Worst Drawdown                  0.566122 0.5635336  0.6442294                 0.601014
Calmar Ratio                    1.437989 1.5831686  0.9208586                 1.367510

My reaction? The execute on next day’s close performance being vastly lower than the other configurations (and that deterioration occurring in the most recent years) essentially means that the fills will have to come pretty quickly at the beginning of the day. While the strategy seems somewhat scalable through the lens of this finger-in-the-air technique, in my opinion, if the first full day of possible execution after signal reception will tank a strategy from a 1.44 Calmar to a .92, that’s a massive drop-off, after holding everything else constant. In my opinion, I think this is quite a valid question to ask anyone who simply sells signals, as opposed to manages assets. Namely, how sensitive are the signals to execution on the next day? After all, unless those signals come at 3:55 PM, one is most likely going to be getting filled the next day.

Now, while this strategy is a bit of a tomato can in terms of how good volatility trading strategies can get (they can get a *lot* better in my opinion), I think it made for a simple little demonstration of this technique. Again, a huge thank you to Mr. Helmuth Vollmeier for so kindly keeping up his dropbox all this time for the volatility data!

Thanks for reading.

NOTE: I am currently contracting in a data science capacity in Chicago. You can email me at, or find me on my LinkedIn here. I’m always open to beers after work if you’re in the Chicago area.

NOTE 2: Today, on October 21, 2015, if you’re in Chicago, there’s a Chicago R Users Group conference at Jaks Tap at 6:00 PM. Free pizza, networking, and R, hosted by Paul Teetor, who’s a finance guy. Hope to see you there.

Volatility Stat-Arb Shenanigans

This post deals with an impossible-to-implement statistical arbitrage strategy using VXX and XIV. The strategy is simple: if the average daily return of VXX and XIV was positive, short both of them at the close. This strategy makes two assumptions of varying dubiousness: that one can “observe the close and act on the close”, and that one can short VXX and XIV.

So, recently, I decided to play around with everyone’s two favorite instruments on this blog–VXX and XIV, with the idea that “hey, these two instruments are diametrically opposed, so shouldn’t there be a stat-arb trade here?”

So, in order to do a lick-finger-in-the-air visualization, I implemented Mike Harris’s momersion indicator.

momersion <- function(R, n, returnLag = 1) {
  momentum <- sign(R * lag(R, returnLag))
  momentum[momentum < 0] <- 0
  momersion <- runSum(momentum, n = n)/n * 100
  colnames(momersion) <- "momersion"

And then I ran the spread through it.

xiv <- xts(read.zoo("longXIV.txt", format="%Y-%m-%d", sep=",", header=TRUE))
vxx <- xts(read.zoo("longVXX.txt", format="%Y-%m-%d", sep=",", header=TRUE))

xivRets <- Return.calculate(Cl(xiv))
vxxRets <- Return.calculate(Cl(vxx))

volSpread <- xivRets + vxxRets
volSpreadMomersion <- momersion(volSpread, n = 252)

In other words, this spread is certainly mean-reverting at just about all times.

And here is the code for the results from 2011 onward, from when the XIV and VXX actually started trading.

#both sides
sig <- -lag(sign(volSpread))
longShort <- sig * volSpread
charts.PerformanceSummary(longShort['2011::'], main = 'long and short spread')

#long spread only
sig <- -lag(sign(volSpread))
sig[sig < 0] <- 0
longOnly <- sig * volSpread
charts.PerformanceSummary(longOnly['2011::'], main = 'long spread only')

#short spread only
sig <- -lag(sign(volSpread))
sig[sig > 0] <- 0
shortOnly <- sig * volSpread
charts.PerformanceSummary(shortOnly['2011::'], main = 'short spread only')

threeStrats <- na.omit(cbind(longShort, longOnly, shortOnly))["2011::"]
colnames(threeStrats) <- c("LongShort", "Long", "Short")
rbind(table.AnnualizedReturns(threeStrats), CalmarRatio(threeStrats))

Here are the equity curves:




With the following statistics:

                          LongShort      Long    Short
Annualized Return          0.115400 0.0015000 0.113600
Annualized Std Dev         0.049800 0.0412000 0.027900
Annualized Sharpe (Rf=0%)  2.317400 0.0374000 4.072100
Calmar Ratio               1.700522 0.0166862 7.430481

In other words, the short side is absolutely amazing as a trade–except for the one small fact of having it be impossible to actually execute, or at least as far as I’m aware. Anyhow, this was simply a for-fun post, but hopefully it served some purpose.

Thanks for reading.

NOTE: I am currently contracting and am looking to network in the Chicago area. You can find my LinkedIn here.

Introduction to Hypothesis Driven Development — Overview of a Simple Strategy and Indicator Hypotheses

This post will begin to apply a hypothesis-driven development framework (that is, the framework written by Brian Peterson on how to do strategy construction correctly, found here) to a strategy I’ve come across on SeekingAlpha. Namely, Cliff Smith posted about a conservative bond rotation strategy, which makes use of short-term treasuries, long-term treasuries, convertibles, emerging market debt, and high-yield corporate debt–that is, SHY, TLT, CWB, PCY, and JNK. What this post will do is try to put a more formal framework on whether or not this strategy is a valid one to begin with.

One note: For the sake of balancing succinctness for blog consumption and to demonstrate the computational techniques more quickly, I’ll be glossing over background research write-ups for this post/strategy, since it’s yet another take on time-series/cross-sectional momentum, except pared down to something more implementable for individual investors, as opposed to something that requires a massive collection of different instruments for massive, institutional-class portfolios.

Introduction, Overview, Objectives, Constraints, Assumptions, and Hypotheses to be Tested:

Momentum. It has been documented many times. For the sake of brevity, I’ll let readers follow the links if they’re so inclined, but among them are Jegadeesh and Titman’s seminal 1993 paper, Mark Carhart’s 1997 paper, Andreu et. Al (2012), Barroso and Santa-Clara (2013), Ilmanen’s Expected Returns (which covers momentum), and others. This list, of course, is far from exhaustive, but the point stands. Formation periods of several months (up to a year) should predict returns moving forward on some holding period, be it several months, or as is more commonly seen, one month.

Furthermore, momentum applies in two varieties–cross sectional, and time-series. Cross-sectional momentum asserts that assets that outperformed among a group will continue to outperform, while time-series momentum asserts that assets that have risen in price during a formation period will continue to do so for the short-term future.

Cliff Smith’s strategy depends on the latter, effectively, among a group of five bond ETFs. I am not certain of the objective of the strategy (he didn’t mention it), as PCY, JNK, and CWB, while they may be fixed-income in name, possess volatility on the order of equities. I suppose one possible “default” objective would be to achieve an outperforming total return against an equal-weighted benchmark, both rebalanced monthly.

The constraints are that one would need a sufficient amount of capital such that fixed transaction costs are negligible, since the strategy is a single-instrument rotation type, meaning that each month may have two-way turnover of 200% (sell one ETF, buy another). On the other hand, one would assume that the amount of capital deployed is small enough such that execution costs of trading do not materially impact the performance of the strategy. That is to say, moving multiple billions from one of these ETFs to the other is a non-starter. As all returns are computed close-to-close for the sake of simplicity, this creates the implicit assumption that the market impact and execution costs are very small compared to overall returns.

There are two overarching hypotheses to be tested in order to validate the efficacy of this strategy:

1) Time-series momentum: while it has been documented for equities and even industry/country ETFs, it may not have been formally done so yet for fixed-income ETFs, and their corresponding mutual funds. In order to validate this strategy, it should be investigated if the particular instruments it selects adhere to the same phenomena.

2) Cross-sectional momentum: again, while this has been heavily demonstrated in the past with regards to equities, ETFs are fairly new, and of the five mutual funds Cliff Smith selected, the latest one only has data going back to 1997, thus allowing less sophisticated investors to easily access diversified fixed income markets a relatively new innovation.

Essentially, both of these can be tested over a range of parameters (1-24 months).

Another note: with hypothesis-driven strategy development, the backtest is to be *nothing more than a confirmation of all the hypotheses up to that point*. That is, re-optimizing on the backtest itself means overfitting. Any proposed change to a strategy should be done in the form of tested hypotheses, as opposed to running a bunch of backtests and selecting the best trials. Taken another way, this means that every single proposed element of a strategy needs to have some form of strong hypothesis accompanying it, in order to be justified.

So, here are the two hypotheses I tested on the corresponding mutual funds:

symbols <- c("CNSAX", "FAHDX", "VUSTX", "VFISX", "PREMX")
getSymbols(symbols, from='1900-01-01')
prices <- list()
for(symbol in symbols) {
  prices[[symbol]] <- Ad(get(symbol))
prices <-, prices)
colnames(prices) <- substr(colnames(prices), 1, 5)
returns <- na.omit(Return.calculate(prices))

sample <- returns['1997-08/2009-03']
monthRets <- apply.monthly(sample, Return.cumulative)

returnRegression <- function(returns, nMonths) {
  nMonthAverage <- apply(returns, 2, runSum, n = nMonths)
  nMonthAverage <- xts(nMonthAverage, = index(returns))
  nMonthAverage <- na.omit(lag(nMonthAverage))
  returns <- returns[index(nMonthAverage)]
  rankAvg <- t(apply(nMonthAverage, 1, rank))
  rankReturn <- t(apply(returns, 1, rank))
  meltedAverage <- melt(data.frame(nMonthAverage))
  meltedReturns <- melt(data.frame(returns))
  meltedRankAvg <- melt(data.frame(rankAvg))
  meltedRankReturn <- melt(data.frame(rankReturn))
  lmfit <- lm(meltedReturns$value ~ meltedAverage$value - 1)
  rankLmfit <- lm(meltedRankReturn$value ~ meltedRankAvg$value)
  return(rbind(summary(lmfit)$coefficients, summary(rankLmfit)$coefficients))

pvals <- list()
estimates <- list()
rankPs <- list()
rankEstimates <- list()
for(i in 1:24) {
  tmp <- returnRegression(monthRets, nMonths=i)
  pvals[[i]] <- tmp[1,4]
  estimates[[i]] <- tmp[1,1]
  rankPs[[i]] <- tmp[2,4]
  rankEstimates[[i]] <- tmp[2,1]
pvals <-, pvals)
estimates <-, estimates)
rankPs <-, rankPs)
rankEstimates <-, rankEstimates)

Essentially, in this case, I take a pooled regression (that is, take the five instruments and pool them together into one giant vector), and regress the cumulative sum of monthly returns against the next month’s return. Also, I do the same thing as the above, except also using cross-sectional ranks for each month, and performing a rank-rank regression. The sample I used was the five mutual funds (CNSAX, FAHDX, VUSTX, VFISX, and PREMX) since their inception to March 2009, since the data for the final ETF begins in April of 2009, so I set aside the ETF data for out-of-sample backtesting.

Here are the results:

pvals <- list()
estimates <- list()
rankPs <- list()
rankEstimates <- list()
for(i in 1:24) {
  tmp <- returnRegression(monthRets, nMonths=i)
  pvals[[i]] <- tmp[1,4]
  estimates[[i]] <- tmp[1,1]
  rankPs[[i]] <- tmp[2,4]
  rankEstimates[[i]] <- tmp[2,1]
pvals <-, pvals)
estimates <-, estimates)
rankPs <-, rankPs)
rankEstimates <-, rankEstimates)

plot(estimates, type='h', xlab = 'Months regressed on', ylab='momentum coefficient', 
     main='future returns regressed on past momentum')
plot(pvals, type='h', xlab='Months regressed on', ylab='p-value', main='momentum significance')
abline(h=.05, col='green')
abline(h=.1, col='red')

plot(rankEstimates, type='h', xlab='Months regressed on', ylab="Rank coefficient",
     main='future return ranks regressed on past momentum ranks', ylim=c(0,3))
plot(rankPs, type='h', xlab='Months regressed on', ylab='P-values')

Of interest to note is that while much of the momentum literature specifies a reversion effect on time-series momentum at 12 months or greater, all the regression coefficients in this case (even up to 24 months!) proved to be positive, with the very long-term coefficients possessing more statistical significance than the short-term ones. Nevertheless, Cliff Smith’s chosen parameters (the two and four month settings) possess statistical significance at least at the 10% level. However, if one were to be highly conservative in terms of rejecting strategies, that in and of itself may be reason enough to reject this strategy right here.

However, the rank-rank regression (that is, regressing the future month’s cross-sectional rank on the past n month sum cross sectional rank) proved to be statistically significant beyond any doubt, with all p-values being effectively zero. In short, there is extremely strong evidence for cross-sectional momentum among these five assets, which extends out to at least two years. Furthermore, since SHY or VFISX, aka the short-term treasury fund, is among the assets chosen, since it’s a proxy for the risk-free rate, by including it among the cross-sectional rankings, the cross-sectional rankings also implicitly state that in order to be invested into (as this strategy is a top-1 asset rotation strategy), it must outperform the risk-free asset, otherwise, by process of elimination, the strategy will invest into the risk-free asset itself.

In upcoming posts, I’ll look into testing hypotheses on signals and rules.

Lastly, Volatility Made Simple has just released a blog post on the performance of volatility-based strategies for the month of August. Given the massive volatility spike, the dispersion in performance of strategies is quite interesting. I’m happy that in terms of YTD returns, the modified version of my strategy is among the top 10 for the year.

Thanks for reading.

NOTE: while I am currently consulting, I am always open to networking, meeting up (Philadelphia and New York City both work), consulting arrangements, and job discussions. Contact me through my email at, or through my LinkedIn, found here.

PELTing a Competing Changepoint Algorithm

This post will demonstrate the PELT algorithm from the changepoint package–a competing algorithm to the twitter package’s breakout detection algorithm. While neither of these algorithms produce satisfactory results, one change point location approximation algorithm that makes no distributional assumptions shows potentially promising results.

I received some feedback regarding my first foray into change point analysis from Brian Peterson. While some of it was good, a fair bit more was how I can add more to my analysis, more boxes I could check off, and so on. One of those boxes was the PELT algorithm, devised by one Rebecca Killick of Lancaster University, which I’ll give a quick run-through of below.

In the twitter paper, PELT was the competing algorithm that the paper compared itself to, and while I didn’t think that replicating the competing algorithm would be necessary at first go, it turns out, that, well, it was necessary. So, going forward, I’m going to have more points demonstrating some more aspects of these change point detection algorithms. Thus far, the most impressive one has been David Matteson’s e.divisive algorithm in his ecp package. However, its one caveat for me is its massively long running time.

Anyhow, without further ado, let’s replicate a diagram found on page 7 of the original paper in the lower right hand corner. Turns out, that is the Scribe data set that comes with the BreakoutDetection package, so we have all we need.

bd <- breakout(Scribe)
pelt <- cpt.meanvar(Scribe)
plot(Scribe, type="l")
abline(v=bd$loc, col="red")
abline(v=pelt@cpts[1], col="blue")
legend(legend = c("Breakout", "PELT"), x = 2, y = 1000, fill = c("red", "blue"))

This gives us the following diagram.

In short, the paper’s claim is corroborated. PELT underperforms even in a simple example, using both packages’ only-one-changepoint methodology. Furthermore, PELT is actually an S4-class type of object (so for those wondering what the @ character is doing, it’s the equivalent of the $ elsewhere in R).

Let’s move onto the GSPC data.


suppressMessages(getSymbols("^GSPC", from = "1984-12-25", to = "2013-05-31"))
#these two lines are only needed if Yahoo gives you holidays such as New Year's -- EG 1985-01-01
GSPC <- GSPC[!as.character(index(GSPC)) %in% as.character(holidayNYSE(1985:2013)),]

dailySqRets <- Return.calculate(Cl(GSPC))*Return.calculate(Cl(GSPC))
dailySqRets <- dailySqRets["1985::"]

Now we’ll look at a histogram of the daily squared returns and see if it looks like some sort of famous distribution.

plot(hist(dailySqRets, breaks=1000), xlim=c(0, 0.001))

Which results in the following image

So, an exponential distribution (give or take)? Well, let’s try and do some PELT changepoint analysis using the exponential distribution.

PELTcps <- cpt.meanvar(as.numeric(dailySqRets), 
                       method="PELT", test.stat="Exponential")

Did that work? No. Here’s the error message.

Error in PELT.meanvar.exp(coredata(data), pen.value) : 
  Exponential test statistic requires positive data

Now, squared returns can’t possibly be negative, because that’s just nonsensical. So what does that mean? Let’s take a look.


And the output:

> dailySqRets[dailySqRets==0]
1985-03-28          0
1985-10-08          0
1988-02-04          0
1988-11-02          0
1992-09-03          0
1997-01-28          0
2003-01-10          0
2008-01-03          0

So, this essentially alleges that there were some days on which the close to close didn’t move at all. Let’s take a look.


This gives us the following output:

> GSPC["1985-03-27::1985-04-01"]
           GSPC.Open GSPC.High GSPC.Low GSPC.Close GSPC.Volume GSPC.Adjusted
1985-03-27    178.43    179.80   178.43     179.54   101000000        179.54
1985-03-28    179.54    180.60   179.43     179.54    99780000        179.54
1985-03-29    179.54    180.66   179.54     180.66   101400000        180.66
1985-04-01    180.66    181.27   180.43     181.27    89900000        181.27

Notice that the close price for the 27th and 28th day are identical, creating a return of zero, which breaks the PELT algorithm. So let’s fix that.

#the method will throw an error with zero returns, so this deals with that
dailySqRets[dailySqRets == 0] <- dailySqRets[dailySqRets==0] + 1e-100 

Essentially, this is just to get past the error messages within the changepoint package. So now, let’s try applying the algorithm once again.

peltCps <- cpt.meanvar(as.numeric(dailySqRets), 
                       method = "PELT", test.stat = "Exponential")

This time, success–it ran! Let’s check the amount of changepoints.


Which gives:

[1] 374

…which is vastly different from the e.divisive algorithm’s output from the previous investigation, or Dr. Robert Frey’s post (20, give or take). The upside is that this algorithm works quickly, but that’s not much solace if the answers are unreliable. To further dive down into the nature of the changepoints, we will remove the last change point (which is simply the length of the data) and do some summary statistics on how long some of these complete “regimes” are (that is, drop the first and last).

newCpts <- peltCps@cpts[-length(peltCps@cpts)]
regimeTime <- diff(newCpts)
hist(regimeTime, breaks = 100)

…which provides the following output:

> summary(regimeTime)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   2.00    3.00    8.00   19.14   26.00  149.00 

And the following histogram:

In short, most of these “regimes” last less than a month, and half of them don’t even make it out to two weeks. These results are not corroborated by the previously investigated methods. As more academic literature uses differences of log returns, and the point is to search for changes in the variance regime, that is the procedure that will be employed, and as the data is continuous and contains negative values, only the Normal distribution is available to choose from when using the PELT method.

logDiffs <- as.numeric(diff(log(Cl(GSPC)))["1985::"])
peltGaussCps <- cpt.var(logDiffs, method="PELT", test.stat="Normal")
fullGaussRegimes <- diff(peltGaussCps@cpts[-length(peltGaussCps@cpts)])
hist(fullGaussRegimes, breaks = 100)

Which gives the following output:

> length(peltGaussCps@cpts)
[1] 93
> fullGaussRegimes <- diff(peltGaussCps@cpts[-length(peltGaussCps@cpts)])
> summary(fullGaussRegimes)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   2.00   16.00   48.00   73.76  106.50  498.00 

And the following histogram:

In short, even though the result is improved, it’s still far from reliable, with most regime switches taking place within 2 months, and many of those lasting barely a few days.

Lastly, I’d like to correct an issue from the previous changepoint post in which I used the “at most one changepoint” method from the BreakoutDetection package. Now, I’ll use the multiple change point detection method.

bdcp <- breakout(dailySqRets, min.size = 20, method = "multi")

With the following results:

> length(bdcp$loc)
[1] 0

In short, the breakout algorithm found no change points whatsoever on default settings, which, once again, does not line up with the results from both the ecp package, along with Dr. Frey’s post. Even if the beta (penalty parameter) gets decreased to .0001 (from .008, its default), it still fails to find any changes in the squared returns data, which is disappointing.

However, over the course of exploring the changepoint package, I have found that there is a method that is directly analogous to Dr. Frey’s cumulative sum of squares method (that is, if you check the help file for cpt.var, one of the test.stat distributions is “CSS”, aka cumulative sum of squares). There are two methods that employ this, neither of which are PELT, but both of which predate PELT (the binary segmentation and segment neighborhood methods), and are explained in Dr. Killick’s original paper.

First off, both algorithms contain a penalty term, and the penalty term that is the default setting is the Bayesian Information Criterion (or BIC aka SIC), which is a double log of the number of data points. In contrast, the AIC is simply the log of 2, so the BIC will be greater than AIC at 8 data points or higher. The cpt.var algorithm function is mostly a wrapper for further nested wrappers, and essentially, the “cut to the chase” is that we eventually nest down to the not-exported binary segmentation variance cumulative sum of squares algorithm, and the segmentation neighborhood sum of squares algorithm.

So here’s the code for the former:

changepoint:::binseg.var.css <- function (data, Q = 5, pen = 0) 
  n = length(data)
  if (n < 4) {
    stop("Data must have atleast 4 observations to fit a changepoint model.")
  if (Q > ((n/2) + 1)) {
    stop(paste("Q is larger than the maximum number of segments", 
               (n/2) + 1))
  y2 = c(0, cumsum(data^2))
  tau = c(0, n)
  cpt = matrix(0, nrow = 2, ncol = Q)
  oldmax = Inf
  for (q in 1:Q) {
    lambda = rep(0, n - 1)
    i = 1
    st = tau[1] + 1
    end = tau[2]
    for (j in 1:(n - 1)) {
      if (j == end) {
        st = end + 1
        i = i + 1
        end = tau[i + 1]
      else {
        lambda[j] = sqrt((end - st + 1)/2) * 
                         ((y2[j + 1] - y2[st])/(y2[end + 1] - y2[st]) - 
                          (j - st + 1)/(end - st + 1))
    k = which.max(abs(lambda))
    cpt[1, q] = k
    cpt[2, q] = min(oldmax, max(abs(lambda), na.rm = T))
    oldmax = min(oldmax, max(abs(lambda), na.rm = T))
    tau = sort(c(tau, k))
  op.cps = NULL
  p = 1:(Q - 1)
  for (i in 1:length(pen)) {
    criterion = (cpt[2, ]) >= pen[i]
    if (sum(criterion) == 0) {
      op.cps = 0
    else {
      op.cps = c(op.cps, max(which((criterion) == TRUE)))
  if (op.cps == Q) {
    warning("The number of changepoints identified is Q, 
             it is advised to increase Q to make sure 
             changepoints have not been missed.")
  return(list(cps = cpt, op.cpts = op.cps, pen = pen))

Essentially, the algorithm loops through all the data points (that is, it defines the start as the beginning of the data, and end as the end of the data), and appends them by the quantity as defined by the lambda[j] line, which is, again:

lambda[j] = sqrt((end - st + 1)/2) * 
                 ((y2[j + 1] - y2[st])/(y2[end + 1] - y2[st]) - 
                  (j - st + 1)/(end - st + 1))

Which, to put it in financial terms, multiplies the square root of half the range by the difference of the percent B (think Bollinger Bands) of the *values* of the data and the percent B of the *indices* of the data for a given start and end location (which are consecutive change point locations). Then, the candidate change point is defined by the maximum of the absolute value of lambda between a starting and ending point, as defined by this code:

    k = which.max(abs(lambda))
    cpt[1, q] = k
    cpt[2, q] = min(oldmax, max(abs(lambda), na.rm = T))
    oldmax = min(oldmax, max(abs(lambda), na.rm = T))
    tau = sort(c(tau, k))

The algorithm then updates tau (the collection of change point locations), which updates the start and end point segment locations, and the algorithm restarts again.

Lastly, at the end, the algorithm loops through the different penalties (in this case, the BIC is simply one constant value, so there may be a special case that allows for some dynamic sort of penalty), and all the change points that have a higher lambda value than the penalty are returned as candidate change points. Again, there is no test of significance in the binary segmentation variance cumulative sum of squares algorithm — so if the penalty is manually specified to be zero, and the user specifies the maximum number of change points (n/2, where n is the length of data), then the algorithm will indeed spit back that many change points. In short, while the default settings may do a half-decent job with finding change points, it’s possible to deliberately force this algorithm to produce nonsensical output. In other words, to be glib, this algorithm isn’t attempting to win the battle against the universe in the never-ending battle to sufficiently idiot-proof something vs. the universe’s capability to create a better idiot. However, using the out-of-the-box settings should sufficiently protect oneself from having absurd results.

Here’s an illustrative example of a few iterations of a few iterations of this algorithm.

In this case, our data will be the daily returns of the GSPC from 1985 onward.

I’ll set Q (the maximum number of change points) at 30.

Let’s start this off.

data <- as.numeric(Return.calculate(Cl(GSPC))["1985::"])
Q <- 30
n <- length(data)
y2 = c(0, cumsum(data^2))
tau = c(0, n)
cpt = matrix(0, nrow = 2, ncol = Q)
oldmax = Inf

Which gives us:

> tau
[1]    0 7328
> cpt
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] 
[1,]    0    0    0    0    0    0    0    0    0    
[2,]    0    0    0    0    0    0    0    0    0  
      [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19]   
[1,]      0     0     0     0     0     0     0     0     0     0
[2,]      0     0     0     0     0     0     0     0     0     0
     [,20] [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30]
[1,]     0     0     0     0     0     0     0     0     0     0     0
[2,]     0     0     0     0     0     0     0     0     0     0     0

In short, tau is the location for the changepoints, and and the cpt matrix will be demonstrated shortly.

Here’s the first iteration through the main loop.

q <- 1
lambda = rep(0, n - 1)
    i = 1
    st = tau[1] + 1
    end = tau[2]
    for (j in 1:(n - 1)) {
      if (j == end) {
        st = end + 1
        i = i + 1
        end = tau[i + 1]
      else {
        lambda[j] = sqrt((end - st + 1)/2) * 
                         ((y2[j + 1] - y2[st])/(y2[end + 1] - y2[st]) - 
                          (j - st + 1)/(end - st + 1))
    k = which.max(abs(lambda))
    cpt[1, q] = k
    cpt[2, q] = min(oldmax, max(abs(lambda), na.rm = T))
    oldmax = min(oldmax, max(abs(lambda), na.rm = T))
    tau = sort(c(tau, k))

And after this first iteration of the loop completes, here are the updated values of tau and cpt:

> tau
[1]    0 5850 7328
> cpt
           [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] 
[1,] 5850.00000    0    0    0    0    0    0    0    0    
[2,]   10.51362    0    0    0    0    0    0    0    0  
      [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19]   
[1,]      0     0     0     0     0     0     0     0     0     0
[2,]      0     0     0     0     0     0     0     0     0     0
     [,20] [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28] [,29] [,30]
[1,]     0     0     0     0     0     0     0     0     0     0     0
[2,]     0     0     0     0     0     0     0     0     0     0     0

In short, tau is the location of the new change point (also the first row of the cpt matrix), and the second row is the absolute value of lambda. After this, the start and end vectors get appended with a new change point to allow for the binary segmentation that forms the basis of the algorithm. However, a picture’s worth a thousand words, so here are some illustrations. The blue lines denote previous change points, and the red point the new change point.

Here is the code I added in the bottom of the loop for plotting purposes:

    k = which.max(abs(lambda))
    plot(lambda, type = "l", main = paste("Changepoint", q))
    abline(v = tau, col="blue")
    cpt[1, q] = k
    cpt[2, q] = min(oldmax, max(abs(lambda), na.rm = T))
    oldmax = min(oldmax, max(abs(lambda), na.rm = T))
    tau = sort(c(tau, k))
    points(x = k, y = lambda[k], col="red", pch = 19)

And here is the result for the first change point:

The second change point:

The 10th change point:

And the 19th (of 30) change points.

Notice that A) lambda is constantly recomputed after every iteration of the main loop, as it’s updated with every new change point and that B) the values of lambda generally decrease as more change point candidates are found, such that the 19th change point is already on the border of the penalty boundary formed by the BIC. Unlike the math in the paper, this particular algorithm in R does not actually stop when lambda of k (that is, lambda[k]) is smaller than the penalty parameter, so if someone plugged in say, Q = 300, the algorithm would take around 30 seconds to run.

So, what’s the punch line to this approximate algorithm?


 t1 <- Sys.time()
 binSegCss <- cpt.var(as.numeric(Return.calculate(Cl(GSPC))["1985::"]), 
                      method="BinSeg", test.stat="CSS", Q = 30)
 t2 <- Sys.time()
 print(t2 - t1)

The algorithm ran in a few seconds for me. Here is the output:

> binSegCss@cpts
 [1]  310  720  739  858 1825 2875 3192 4524 4762 5044 5850 6139 6199 6318 6548 6641 6868 6967 7328
> length(binSegCss@cpts)
[1] 19

Since the last change point is the length of the data, we’ll disregard that. In short, 18 change points (so that last picture of the four change points? That one and all subsequent ones fell in the realm of noise), which falls right into the ballpark of the post from Keplerian Finance.

So, as before, let’s create the cumulative sum of squares plot, the time series plot, and the volatility terrain map.

returns <- Return.calculate(Cl(GSPC))["1985::"]
dailySqReturns <- returns*returns
cumSqReturns <- cumsum(dailySqReturns)
xtsExtra::addLines(index(dailySqReturns)[binSegCss@cpts[-length(binSegCss@cpts)]], on = 1, col="blue", lwd = 2)

xtsExtra::addLines(index(dailySqReturns)[binSegCss@cpts[-length(binSegCss@cpts)]], on = 1, col="blue", lwd = 2)

returns$regime <- NA
for(i in 1:(length(binSegCss@cpts)-1)) {
  returns$regime[binSegCss@cpts[i]] <- i
returns$regime <- na.locf(returns$regime)
returns$regime[$regime)] <- 0
returns$regime <- returns$regime + 1
returns$annVol <- NA
for(i in 1:max(unique(returns$regime))) {
  regime <- returns[returns$regime==i,]
  annVol <- StdDev.annualized(regime[,1])
  returns$annVol[returns$regime==i,] <- annVol


Which gives us the following three images:

This last volatility map is even closer to the one found in Keplerian Finance, thus lending credibility to the technique.

In conclusion, it seems that the twitter breakout detection algorithm (e-divisive with medians) is “too robust” against the type of events in financial data, and thus does not detect enough change points. On the other hand, the PELT algorithm suffers from the opposite issues — by forcing the assumption of a popular distribution from probability theory (EG Normal, Exponential, Poisson, Gamma), PELT outputs far too many candidate change points, making its results highly suspect. However, a special case of the binary segmentation algorithm — the binary segmentation algorithm for variance using cumulative sum of squares — presents interesting results. Although the algorithm is a heuristic approximation, its fast running time and lack of distributional assumptions lend it usefulness for doing a “quick and dirty” analysis to potentially cut down on the running time of a more computationally-intensive changepoint detecting algorithm.

When this topic will be revisited, the segment neighborhood method will be examined.

Thanks for reading.

NOTE: I am a freelance consultant in quantitative analysis on topics related to this blog. If you have contract or full time roles available for proprietary research that could benefit from my skills, please contact me through my LinkedIn here.