CorrelationFunction for vectors

Question

Mathematica 9 has a new CorrelationFunction. Sadly the joy does not last long, as I can't get it to work with vectors. I would like to make a velocity correlation function. That is, given a list of n vectors I would like to calculate:

For example, this is almost instantaneous:

testdata = RandomVariate[BinormalDistribution[{10, 20}, 0.4], {10000}];
AbsoluteTiming[
x = CorrelationFunction[testdata[[All, 1]], {0, 1000}];
y = CorrelationFunction[testdata[[All, 2]], {0, 1000}];
]
ListPlot[{x, y}, PlotRange -> {0, 1}]

^{(the test data here is not the best as it is by definition uncorrelated)}

Here I would like the autocorrelation of 2D vectors evenly spaced in time. But this does not work:

r = CorrelationFunction[testdata, {0, 100}];

How can I make an autocorrelation for vectors? Or do I have to go the manual way?

ClearAll[AbsoluteAutocorrelationT, AutocorrelationT]

AbsoluteAutocorrelationT[l_, 0] := Total[(#.#) & /@ l];
AbsoluteAutocorrelationT[l_, n_] := 
  Total@Table[l[[i]].l[[i + n]], {i, Length[l] - n}];
AutocorrelationT[l_, range_] := Block[{t0, ti},
   t0 = AbsoluteAutocorrelationT[l, 0];
   ti = AbsoluteAutocorrelationT[l, #] & /@ (Range @@ range);
   (ti/t0)];

AbsoluteTiming[r = AutocorrelationT[testdata, {1, 1000}];]

Which is 200 times slower. I know I can compile this (but than I have to worry about overflows in the sum, etc.)

Answer 1

When each $v_i = (v_{i,1}, v_{i,2}, \ldots, v_{i,k})$ is a $k$-vector we may interchange the summations to obtain

$$\sum_{i=1}^{n-t} v_i\cdot v_{i+t} = \sum_{j=1}^k\sum_{i=1}^{n-t} v_{i,j} v_{i+t,j},$$

whence the desired correlation is the sum of the correlations of the components. Thus:

kernel = PDF[NormalDistribution[0, 1]] /@ Range[-5, 5, .005]; (* Induces autocorrelation *)
testdata = ListConvolve[kernel, #] & /@
  Transpose[RandomVariate[BinormalDistribution[{10, 20}, 0.4], {10^4}]];
c = Sum[ListCorrelate[v, v, {1, 1}], {v, testdata}]; // AbsoluteTiming

$\{0.0030002, \text{Null}\}$

ListLinePlot[c[[;; Length[c]/2]]/c[[1]]]

The timings scale as expected ($k=2$ shown):

Answer 2

TL/DR:

temporaltestData = TemporalData[{testdata},Automatic];
r = CorrelationFunction[temporaltestData, {0, 1000}]["Values"];

---

Recap of Question

In your question, you have provided the code

testdata = RandomVariate[BinormalDistribution[{10, 20}, 0.4], {10000}];
AbsoluteTiming[
  x = CorrelationFunction[testdata[[All, 1]], {0, 1000}];
  y = CorrelationFunction[testdata[[All, 2]], {0, 1000}];
]
ListPlot[{x, y}, PlotRange -> {0, 1}]

and you found that this code finds the autocorrelations in the first and second components of your time series just fine. However, you also wanted to find the lagged cross-correlations between the first and second components.

Attempted solution

You attempted to do this with

r = CorrelationFunction[testdata, {0, 100}];

but you found that this doesn't work. I got an error that said

CorrelationFunction::tmpln: Uniformly spaced sampling expected in {{-3.38558,-24.5187},{-5.09168,-0.21692},{7.14122,14.2209},{-9.21681,19.2432},{8.15215,14.9583},{0.911419,-8.53165},{-19.1743,-22.1},{10.885,14.1329},{5.9896,3.49815},{9.51289,32.5084},<<31>>,{-1.15955,25.1969},{11.6961,13.6815},{-13.5907,-14.7231},{15.2179,17.0268},{4.61082,11.2426},{4.04496,-2.09691},{-0.278494,17.8067},{-4.02344,10.879},{19.5057,8.53984},<<9950>>}.

What I believe is happening is that mathematica is interpreting the first component as a time. You can try to get around this by prepending the times:

timedTestData = Prepend[testdata[[#]], #] & /@ Range@Length@testdata
r = CorrelationFunction[timedTestData, {0, 100}]

but now I get the error

CorrelationFunction::bdlag: The lag specification {0,100} should be a symbol, an integer with magnitude less than the length of the data, or a range specification indicating such integers.

which doesn't make sense.

Real solution

However, there is a way to find the correlation, as you can see by going to the help for CorrelationFunction, and then going to Examples->Scope->Empirical Estimates, and looking at the last example in this section. The example is

proc = ARProcess[{{{.3, .1}, {.6, .3}}}, {{1, .3}, {.3, .6}}];
data = RandomFunction[proc, {100}];
cov = CorrelationFunction[data, {-6, 6}];

They have no problem taking the correlation of the multidimensional data set. This is because their data is a TemporalData object. So we must make one of our own. This can be done with

temporaltestData = TemporalData[{testdata},Automatic];

The brackets around testdata are necessary because otherwise TemporalData will treat testdata as 10000 time series of length two instead of one 2D time series of length 10000. The second argument, Automatic, just assigns consecutive integer times to each data point. Once the data is packaged in the TemporalData object, CorrelationFunction has a much easier time knowing what to do with it:

r = CorrelationFunction[temporaltestdata, {0, 1000}]["Values"]

Here, CorrelationFunction returns a TimeSeries object, which is a special kind of TemporalData object. If we just want the list of lagged autocovariance matrices, we can access them by applying our time series object to "Values". Now we have our lagged autocovariances. Notice for example:

Norm@Chop[r[[;; , 1, 1]] - x]

gives 0

Answer 3

The following is (only) six times faster:

l = RandomVariate[BinormalDistribution[{10, 20}, 0.4], {10000}];

c[l_, 0] := ListCorrelate[l, l, {1, -1}, {}, Dot, Plus, 1];
c[l_, s_List] := Module[{m = c[l, 0][[1]]}, 
   ListCorrelate[l[[;;Length@l - #]], l[[1+#;;]], {1, -1}, {}, Dot, Plus, 1] / m & /@ Range @@ s];

ListLinePlot@Table[ AbsoluteTiming[AutocorrelationT[l, {1, n}];][[1]]/
                    AbsoluteTiming[c[l, {1, n}];][[1]], {n, 30}]