Ehlers’s Autocorrelation Periodogram

This post will introduce John Ehlers’s Autocorrelation Periodogram mechanism–a mechanism designed to dynamically find a lookback period. That is, the most common parameter optimized in backtests is the lookback period.

Before beginning this post, I must give credit where it’s due, to one Mr. Fabrizio Maccallini, the head of structured derivatives at Nordea Markets in London. You can find the rest of the repository he did for Dr. John Ehlers’s Cycle Analytics for Traders on his github. I am grateful and honored that such intelligent and experienced individuals are helping to bring some of Dr. Ehlers’s methods into R.

The point of the Ehlers Autocorrelation Periodogram is to dynamically set a period between a minimum and a maximum period length. While I leave the exact explanation of the mechanic to Dr. Ehlers’s book, for all practical intents and purposes, in my opinion, the punchline of this method is to attempt to remove a massive source of overfitting from trading system creation–namely specifying a lookback period.

SMA of 50 days? 100 days? 200 days? Well, this algorithm takes that possibility of overfitting out of your hands. Simply, specify an upper and lower bound for your lookback, and it does the rest. How well it does it is a topic of discussion for those well-versed in the methodologies of electrical engineering (I’m not), so feel free to leave comments that discuss how well the algorithm does its job, and feel free to blog about it as well.

In any case, here’s the original algorithm code, courtesy of Mr. Maccallini:

AGC <- function(loCutoff = 10, hiCutoff = 48, slope = 1.5) {      accSlope = -slope # acceptableSlope = 1.5 dB   ratio = 10 ^ (accSlope / 20)   if ((hiCutoff - loCutoff) > 0)
    factor <-  ratio ^ (2 / (hiCutoff - loCutoff));
  return (factor)
}

autocorrPeriodogram <- function(x, period1 = 10, period2 = 48, avgLength = 3) {
  # high pass filter
  alpha1 <- (cos(sqrt(2) * pi / period2) + sin(sqrt(2) * pi / period2) - 1) / cos(sqrt(2) * pi / period2)
  hp <- (1 - alpha1 / 2) ^ 2 * (x - 2 * lag(x) + lag(x, 2))
  hp <- hp[-c(1, 2)]
  hp <- filter(hp, (1 - alpha1), method = "recursive")
  hp <- c(NA, NA, hp)
  hp <- xts(hp, order.by = index(x))
  # super smoother
  a1 <- exp(-sqrt(2) * pi / period1)
  b1 <- 2 * a1 * cos(sqrt(2) * pi / period1)
  c2 <- b1
  c3 <- -a1 * a1
  c1 <- 1 - c2 - c3
  filt <- c1 * (hp + lag(hp)) / 2
  leadNAs <- sum(is.na(filt))
  filt <- filt[-c(1: leadNAs)]
  filt <- filter(filt, c(c2, c3), method = "recursive")
  filt <- c(rep(NA, leadNAs), filt)
  filt <- xts(filt, order.by = index(x))
  # Pearson correlation for each value of lag
  autocorr <- matrix(0, period2, length(filt))
  for (lag in 2: period2) {
    # Set the average length as M
    if (avgLength == 0) M <- lag
    else M <- avgLength
    autocorr[lag, ] <- runCor(filt, lag(filt, lag), M)
  }
  autocorr[is.na(autocorr)] <- 0
  # Discrete Fourier transform
  # Correlate autocorrelation values with the cosine and sine of each period of interest
  # The sum of the squares of each value represents relative power at each period
  cosinePart <- sinePart <- sqSum <- R <- Pwr <- matrix(0, period2, length(filt))
  for (period in period1: period2) {
    for (N in 2: period2) {
      cosinePart[period, ] = cosinePart[period, ] + autocorr[N, ] * cos(2 * N * pi / period)
      sinePart[period, ] = sinePart[period, ] + autocorr[N, ] * sin(2 * N * pi / period)
    }
    sqSum[period, ] = cosinePart[period, ] ^ 2 + sinePart[period, ] ^ 2
    R[period, ] <- EMA(sqSum[period, ] ^ 2, ratio = 0.2)
  }
  R[is.na(R)] <- 0
  # Normalising Power
  K <- AGC(period1, period2, 1.5)
  maxPwr <- rep(0, length(filt))   for(period in period1: period2) {     for (i in 1: length(filt)) {       if (R[period, i] >= maxPwr[i]) maxPwr[i] <- R[period, i]
      else maxPwr[i] <- K * maxPwr[i]
    }
  }
  for(period in 2: period2) {
    Pwr[period, ] <- R[period, ] / maxPwr
  }
  # Compute the dominant cycle using the Center of Gravity of the spectrum
  Spx <- Sp <- rep(0, length(filter))
  for(period in period1: period2) {
    Spx <- Spx + period * Pwr[period, ] * (Pwr[period, ] >= 0.5)
    Sp <- Sp + Pwr[period, ] * (Pwr[period, ] >= 0.5)
  }
  dominantCycle <- Spx / Sp
  dominantCycle[is.nan(dominantCycle)] <- 0
  dominantCycle <- xts(dominantCycle, order.by=index(x))
  dominantCycle <- dominantCycle[dominantCycle > 0]
  return(dominantCycle)
  #heatmap(Pwr, Rowv = NA, Colv = NA, na.rm = TRUE, labCol = "", add.expr = lines(dominantCycle, col = 'blue'))
}

One thing I do notice is that this code uses a loop that says for(i in 1:length(filt)), which is an O(data points) loop, which I view as the plague in R. While I’ve used Rcpp before, it’s been for only the most basic of loops, so this is definitely a place where the algorithm can stand to be improved with Rcpp due to R’s inherent poor looping.

Those interested in the exact logic of the algorithm will, once again, find it in John Ehlers’s Cycle Analytics For Traders book (see link earlier in the post).

Of course, the first thing to do is to test how well the algorithm does what it purports to do, which is to dictate the lookback period of an algorithm.

Let’s run it on some data.

getSymbols('SPY', from = '1990-01-01')

t1 <- Sys.time()
out <- autocorrPeriodogram(Ad(SPY), period1 = 120, period2 = 252, avgLength = 3)
t2 <- Sys.time() print(t2-t1) 

And the result:

 > t1 <- Sys.time() > out <- autocorrPeriodogram(Ad(SPY), period1 = 120, period2 = 252, avgLength = 3) > t2 <- Sys.time() > print(t2-t1)
Time difference of 33.25429 secs

Now, what does the algorithm-set lookback period look like?

plot(out)

Let’s zoom in on 2001 through 2003, when the markets went through some upheaval.

plot(out['2001::2003']

In this zoomed-in image, we can see that the algorithm’s estimates seem fairly jumpy.

Here’s some code to feed the algorithm’s estimates of n into an indicator to compute an indicator with a dynamic lookback period as set by Ehlers’s autocorrelation periodogram.

acpIndicator <- function(x, minPeriod, maxPeriod, indicatorFun = EMA, ...) {
  acpOut <- autocorrPeriodogram(x = x, period1 = minPeriod, period2 = maxPeriod)
  roundedAcpNs <- round(acpOut, 0) # round to the nearest integer
  uniqueVals <- unique(roundedAcpNs) # unique integer values
  out <- xts(rep(NA, length(roundedAcpNs)), order.by=index(roundedAcpNs))

  for(i in 1:length(uniqueVals)) { # loop through unique values, compute indicator
    tmp <- indicatorFun(x, n = uniqueVals[i], ...)
    out[roundedAcpNs==uniqueVals[i]] <- tmp[roundedAcpNs==uniqueVals[i]]
  }
  return(out)
}

And here is the function applied with an SMA, to tune between 120 and 252 days.

ehlersSMA <- acpIndicator(Ad(SPY), 120, 252, indicatorFun = SMA)

plot(Ad(SPY)['2008::2010'])
lines(ehlersSMA['2008::2010'], col = 'red')

And the result:

As seen, this algorithm is less consistent than I would like, at least when it comes to using a simple moving average.

For now, I’m going to leave this code here, and let people experiment with it. I hope that someone will find that this indicator is helpful to them.

Thanks for reading.

NOTES: I am always interested in networking/meet-ups in the northeast (Philadelphia/NYC). Furthermore, if you believe your firm will benefit from my skills, please do not hesitate to reach out to me. My linkedin profile can be found here.

Lastly, I am volunteering to curate the R section for books on quantocracy. If you have a book about R that can apply to finance, be sure to let me know about it, so that I can review it and possibly recommend it. Thakn you.