I have a pretty long time series (4000 data points).

How can I calculate the Correlation dimension and/or Lyapunov exponent for such data?

Please also advise on how to import the data as well.

  • 1
    $\begingroup$ What is the format of the data? Depending on the sort of file it is in you can often simply define datain=Import["/yourdirectory/yourfile.dat"] but you should look at the options for Import which may help. $\endgroup$ Commented May 28, 2013 at 10:29
  • 1
    $\begingroup$ It will be also helpful if you supply us some example data to play with. You can upload your data somewhere. $\endgroup$ Commented May 28, 2013 at 10:30
  • 3
    $\begingroup$ I fail to see down-voter's reasoning. $\endgroup$
    – Sasha
    Commented May 28, 2013 at 13:55

3 Answers 3


(EDIT: Another implementation has been added at the end of this post.)

I employ Nearest to implement the algorithm for correlation dimension, $d_C$. Here's how it goes:

Make some data

data = RandomReal[1, {4000, 2}];
ListPlot[data, Frame -> True, PlotStyle -> Black, AspectRatio -> 1]

enter image description here

I define the correlation sum using Nearest

corr[r_] := Total[Table[
Length[Drop[Nearest[Drop[data, i - 1], data[[i]], {All, r}], 
 1]], {i, 1, Length[data]}]]

The Drop between Length and Nearest is because Nearest[dataset,dataset[[i]]] results in dataset[[i]] for any i.

For a comparison I define also the correlation sum in the regular way

corrSum[r_] := 
Sum[HeavisideTheta[r - EuclideanDistance[data[[i]], data[[j]]]], {i, 
1, Length[data]}, {j, i + 1, Length[data]}]

(I don't care about the multiplicative constant in front of these sums because later I will take a logarithm which will result in just a vertical offset not influencing the slope and hence the correlation dimension.)

I evaluate the two approaches:

my = Log[10, N@Table[{r, corr[r]}, {r, 0.001, 0.01, 0.001}]]; // AbsoluteTiming

{1.83536, Null}

classic = Log[10, N@Table[{r, corrSum[r]}, {r, 0.001, 0.01, 0.001}]]; // AbsoluteTiming

{187.483, Null}

My implementation is 100 times faster than the straightforward approach, and yields the correct result:

my == classic


Now only some plotting and fitting:

plot1 = ListPlot[my, PlotStyle -> {Black, PointSize[Large]}];
fit = LinearModelFit[Drop[my, 1], x, x];
{min, max} = MinMax@First@Transpose@my;
plot2 = Plot[Normal[fit], {x, min, max}, PlotStyle -> Red];
Show[plot1, plot2, Frame -> True, PlotRange -> All, 
PlotLabel -> "corr_dim=" <> ToString@Normal[fit][[2, 1]], 
FrameLabel -> {"log(r)", "log[C(r)]"}]

to achieve

enter image description here

The idea is to fit a line in the linear part of the $\log C(r)$ vs. $\log r$ region, hence for the dataset data I had to Drop[my,1].

The correlation dimension $D_C$ will usually give a value slightly smaller than the correct fractal dimension; here, $d_C=1.96$ instead of 2, which is not that bad.

The set of values of r that I used above is only exemplary and should be adjusted in a particular case. The same applies to chosing a subset of my for fitting a line.

I ran the exact same code on the Duffing attractor (also 4000 points)

enter image description here

and achieved

enter image description here

i.e. $d_C=1.33$; time elapsed

{5.5014, Null}

EDIT: This is directly inspired by an answer given by C. E. to a different question. In fact, I admit I stole even his names of symbols and functions.

This is a SparseArray approach, so I refer to it as SA:

corrSA[data_, r_] := Module[{dist, eliminate, pos},
  dist = DistanceMatrix[data, DistanceFunction -> SquaredEuclideanDistance];
  eliminate = 
   UnitStep[ConstantArray[r^2, {Length@data, Length@data}] - dist] - 
  pos = SparseArray[eliminate]["NonzeroPositions"];
  0.5 Length@pos

Note that I employed DistanceFunction -> SquaredEuclideanDistance in dist because it turns out it is faster than EuclideanDistance. Accordingly, I had to put r^2 in eliminate.

Sadly, for random data this takes more time than corr[r] took:

mySA = Log10[N@Table[{r, corrSA[data, r]}, {r, 0.001, 0.01, 0.001}]]; // AbsoluteTiming

{2.837545, Null}

On the other hand, it required substantially less time for the Duffing attractor:

{2.81612, Null}

EDIT2: All computations were performed without any parallelization. In cases when ParallelTable was used to compute my and classic (i.e., employing the definitions in corr and corrSum, respectively), a 3x speed-up was achieved on a 4-core machine.

In case of obtaining mySA with parallelization, the speed-up was only by 14%.


Appendix 1 of Mario Martelli's Introduction to Discrete Dynamical Systems and Chaos contains sample Mathematica code (pages 290-296) which should help you http://dx.doi.org/10.1002/9781118032879.app1 (free, click 'Get PDF')

Also, the Wolfram Library Archive contains sample code in Testing Chaos and Fractal Properties in Economic Time Series computes these quantities for 1,000 data points. http://library.wolfram.com/infocenter/Conferences/6162/

  • $\begingroup$ In case anyone wants to follow the second lead, the linked file seems inaccessible, but you can get it through the Internet Wayback Machine at (web.archive.org/web/20100821131053/http://…). It's a bit hard to understand, but it looks like they might be implementing the Wolf et al. 1985 algorithm for calculating Lyapunov exponents from a time series. $\endgroup$
    – Chris K
    Commented Aug 26, 2016 at 20:29
  • $\begingroup$ That Wolf et al. 1985 algorithm has also been ported to Matlab. Hope someone follows up on this to provide a nice answer! $\endgroup$
    – Chris K
    Commented Aug 26, 2016 at 20:30

An alternative way of using Nearest gives over 15X speed-up over corrSA in corey979's answer:

corrNearest[data_, r_] := Module[{nF = Nearest[data -> "Index"]},
  Total[Table[UnitStep[nF[data[[i]], {All, r}]- i - 1], {i, 1, Length[data]}], 2]] 


n = 2000;
data = RandomReal[1, {n, 2}];
res = Log10 @ N @ Table[{r, corrNearest[data, r]}, {r, 0.001, 0.01, 0.001}]; //
  AbsoluteTiming // First


mySA = Log10 @ N @ Table[{r, corrSA[data, r]}, {r, 0.001, 0.01, 0.001}]; // 
  AbsoluteTiming // First


res == mySA 


For n = 3000 the timings are 2.75211 for corrSA and 0.132438 for corrNearest.

  • $\begingroup$ Which version of MMA does this apply to? MMA (win10 64 bit) seems not to like data->"Index" (and tbh I don't understand what that is doing). I get "Index is neither a list of real points not a valid list of rules" when I use corrNearest (Code just copy/pasted...) $\endgroup$ Commented Jun 15, 2019 at 10:46
  • $\begingroup$ @Julian, I think it was v10 (wolfram cloud) when this was posted. Just tried teh same code in version 12 (wolfram cloud) and it works. The form data -> "Index" is introduced in version 10; it should work in later versions. The form data ->Automatic works in version 9 and newer versions and it does the same thing, that is, Nearest[data ->"Index"][x] and Nearest[data ->Automatic][x] both give the position(s) of the element in data that is closest to x. $\endgroup$
    – kglr
    Commented Jun 15, 2019 at 10:59
  • $\begingroup$ I see you also said it was available in 11 in this Q mathematica.stackexchange.com/questions/25273/… but it does not seem to work for me. (checked with fresh kernel) However data->Automatic does. Thanks $\endgroup$ Commented Jun 15, 2019 at 11:03
  • $\begingroup$ @Julian, it is probably a glitch introduced in the MMA release fixed later. $\endgroup$
    – kglr
    Commented Jun 15, 2019 at 11:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.