Validate Cross Correlation, Part 2

If you are unfamiliar with the markets, or how to interpret a market feed as a real function you might want to check the previous post on this topic.

Let’s start with a simple function and apply the fourier transform and its inverse, the R snippets are available in the repository if you prefer to see or modify the code. Here we will break after each bit of code to offer some explanations.

first we write a simple function to generate triangular functions. Nothing fancy really, but will save us time later

triangle <- function(period) {
  p4 <- period / 4
  up <- (seq(1, period/2) - p4) / p4
  dn <- (p4 - seq(1, period/2)) / p4
  return(c(up, dn))
}

using the function we create a triangle

t <- triangle(128)

and wrap the triangle in a data.frame(), because ggplot2 really likes data.frame()

df <- data.frame(usec=seq(0, length(t) - 1), value=t, variable="T(t)")

ggplot2 generates sensible and good looking plots in most cases, make sure it is loaded

require(ggplot2)  # May need to install.packages("ggplot2")

then we can plot the triangle function

qplot(x=usec, y=value, color=variable, data=df)

next, let’s apply the FFT transform and the inverse to the triangular function

fft.i <- fft(fft(t), inverse=TRUE)

and save this into a new data.frame()

tmp <- df
tmp$value <- Re(fft.i)
tmp$variable <- 'FFT^1(FFT(t))'
df.i <- rbind(df, tmp)

let’s add the new data to the data.frame()

qplot(x=usec, y=value, color=variable, data=df.i)

What is going on here? To save computations, the “Fast” Fourier Transform omits rescaling the function by $1/N$, where $N$ is the number of samples. If we apply this rescaling manually things match perfectly

fft.i <- fft.i / length(fft.i)
tmp <- df
tmp$value <- Re(fft.i)
tmp$variable <- 'FFT^1(FFT(T))(t)'
df.i <- rbind(df, tmp)
qplot(x=usec, y=value, color=variable, data=df.i)

the more or less obvious question is how well does a function correlate to itself, this is easy to compute

corr <- Re(fft( Conj(fft(t)) * fft(t), inverse=TRUE)) / length(t)

we are going to be wrapping these functions in a data.frame() a lot, so let’s create a function for it

to.df <- function(a, name) {
  return(data.frame(
    usec=seq(0, length(a) - 1), value=a, variable=name))
}

and plot the results

corr.df <- to.df(corr, "(T * T)(t)")
qplot(x=usec, y=value, color=variable, data=corr.df) +
  scale_y_continuous(name=expression(value^2)) +
  theme(legend.position="bottom")

Notice that the y-axis is labeled $value^2$, this is because the result of a cross correlation is a measure of area. That means that as we process data with larger values the cross-correlation will grow with the square of the value too.

Let’s see how the correlation works with a time shifted signal

a <- triangle(128)
a.df <- to.df(a, "A")
b <- a[((seq(1, length(a)) - 13) %% length(a)) + 1]
ab.df <- rbind(a.df, to.df(b, "B"))
qplot(x=usec, y=value, color=variable, data=ab.df) +
  theme(legend.position="bottom")

Let’s see how the cross-correlation looks like, but since we will be doing several correlations, we write a helper function …

correlation <- function(a, b) {
  inv <- Re(fft( Conj(fft(a)) * fft(b), inverse=TRUE))
  return(inv / length(a))
}
correlation.df <- function(a, b, name) {
  return(to.df(correlation(a, b), name))
}
corr.ab.df <- correlation.df(a, b, "A * B")
qplot(x=usec, y=value, color=variable, data=corr.ab.df) +
  scale_y_continuous(name=expression(value^2)) +
  theme(legend.position="bottom")

The graphs are pretty, but exactly where is the peak?

which.max(corr.ab.df$value)

[1] 13

That is a perfect match, but market (or other) signals are rarely so perfectly match, what happens if we add some noise?