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.

Pingback: Ehlers’s Autocorrelation Periodogram – Mubashir Qasim

Pingback: Quantocracy's Daily Wrap for 02/15/2017 | Quantocracy

Pingback: Ehlers’s Autocorrelation Periodogram | Econometria aplicada

Pingback: Distilled News | Data Analytics & R

Excellent as always

This emits multiple errors in R:

maxPwr = maxPwr[i]) maxPwr[i] <- R[period, i]

Error: unexpected 'for' in:

" K <- AGC(period1, period2, 1.5)

maxPwr else maxPwr[i] }

Error: unexpected ‘}’ in ” }”

> }

Error: unexpected ‘}’ in ” }”

> for(period in 2: period2) {

+ Pwr[period, ] # Compute the dominant cycle using the Center of Gravity of the spectrum

> Spx <- Sp for(period in period1: period2) {

+ Spx = 0.5)

+ Sp = 0.5)

+ }

Error: object ‘period1’ not found

> dominantCycle dominantCycle[is.nan(dominantCycle)] dominantCycle dominantCycle 0]

> return(dominantCycle)

Error: no function to return from, jumping to top level

> #heatmap(Pwr, Rowv = NA, Colv = NA, na.rm = TRUE, labCol = “”, add.expr = lines(dominantCycle, col = ‘blue’))

> }

Error: unexpected ‘}’ in “}”

It worked for me. I copied and pasted the code from the github.

This was fine, however

https://github.com/fabmaccallini/ehlers/blob/master/autocorrPeriodogram.R

cant seem to plot (out)

it says

plot(out)

Error in xy.coords(x, y, xlabel, ylabel, log) :

‘x’ is a list, but does not have components ‘x’ and ‘y’

Idea here is that stock prices are random with memory. Hence autocorrelation picks up memory in the series. If you have a perfect sine wave with peak and trough every 10 bars. At the cycle of 10 bars you see non correlation and on 20 bars you see correlation between points. The switching of the correlation coefficients at the full 20 bar cycle period and 10 bar cycle period mark the turning points.