Using mdDelay() and mdFnn() to estimate embedding of the Lorenz attractor

This file contains examples of how to use the functions mdDelay and mdFnn on example data from the Lorenz equations. This script was used to produce Figure 1 and Figure 2 in the article, but also contains some additional examples of calling the functions and plotting the results.

For further details, we refer to the article:

Wallot, S., & Mønster, D. (2018). Calculation of average mutual information (AMI) and false-nearest neighbors (FNN) for the estimation of embedding parameters of multidimensional time-series in Matlab. Front. Psychol. - Quantitative Psychology and Measurement (under review)

Contents

Load the data and set font size for plots

data = load('lorenz_3d_timeseries.txt');
fontSize = 18;

Plot the full Lorenz attractor (Figure 1a in the article)

figure1a = figure;
set(gcf,'color','white')
axes1a = axes('Parent',figure1a);
hold(axes1a,'on');
plot3(data(:,1),data(:,2),data(:,3),'k')
xlabel('x')
ylabel('y')
zlabel('z')
view(axes1a,[67 25]);
set(gca,'FontSize',fontSize,'fontWeight','normal')
print('Figure1a','-dpng')

Estimate the time delay for the x-variable

The estimate returned here is 19, which in this case is higher than it needs to be. This can bee seen from the plot, where the AMI function does not drop below the threshold value of 1/e. The function thus falls back to finding the first local minimum which occurs at a delay of 19. Since the AMI function is very flat (the first derivative is almost zero), a value very close to the first local minimum occurs already at a delay of 13 or 14, so there is no reason to choose a larger value than that.

This shows the importance of visually inspecting the results. For y and z the algorithm produces values that hold up to visual inspection, and the same is true for the combinations xy, xz, yz and xyz (see plots further below).

tau = mdDelay(data(:,1), 'maxLag', 25, 'plottype', 'all');
set(gca,'FontSize',fontSize,'fontWeight','normal')
disp('x: tau = ' + string(tau))

No value below threshold found. Will use first local minimum instead
x: tau = 19

Construct time delayed versions of the x-variable (Figure 1b-d in the article)

To see the effect in the plot we use an exaggerated value of tau, so for illustration purposes we set tau = 40.

tau = 40;
figure();
set(gcf,'color','white')
% Plot the x variable
subplot(3,1,1), plot(data(:,1),'k')
ylabel('x')
axis([0 2000 min(data(:,1)) max(data(:,1))])
set(gca,'FontSize',fontSize,'fontWeight','normal')
set(findall(gcf,'type','text'),'FontSize',fontSize,'fontWeight','normal')
% Plot the x variable delayed by tau
subplot(3,1,2), plot(data(1 + tau:end,1),'k')
ylabel('x')
axis([0 2000 min(data(:,1)) max(data(:,1))])
set(gca,'FontSize',fontSize,'fontWeight','normal')
set(findall(gcf,'type','text'),'FontSize',fontSize,'fontWeight','normal')
subplot(3,1,3), plot(data(1 + 2 * tau:end,1),'k')
ylabel('x')
xlabel('time')
axis([0 2000 min(data(:,1)) max(data(:,1))])
set(gca,'FontSize',fontSize,'fontWeight','normal')
set(findall(gcf,'type','text'),'FontSize',fontSize,'fontWeight','normal')

Reconstruct the attractor based on the x time series (Figure 1e in the article)

The x time series is embedded using time delayed versions of the time series with a time delay tau and embedding dimension m. As discussed above, visual inspection of the AMI function indicates that the first local minimium value is too high a value for tau, since the AMI reaches essentially the same level at tau around 13, so we will pick this value. The embedding dimension is set to 3, which is the value produced by mdFnn() (see Table 1 in the article).

tau = 13;
m = 3;
figure1e = figure();
set(gcf,'color','white')
axes1e = axes('Parent',figure1e);
hold(axes1e,'on');
plot3(data(1:end-(m-1)*tau,1),data(1+tau:end-(m-2)*tau,1),data(1+2*tau:end,1),'k')
xlabel('x')
ylabel('y')
zlabel('z')
view(axes1e,[67 25]);
set(gca,'FontSize',fontSize,'fontWeight','normal')
set(findall(gcf,'type','text'),'FontSize',fontSize,'fontWeight','normal')

Estimate time delay and plot AMI using all variables (Figure 2a in article)

tau = mdDelay(data, 'maxLag', 25, 'plottype', 'all');
set(gca,'FontSize',fontSize,'fontWeight','normal')
disp('xyz: tau = ' + string(tau))
print('Figure2a','-dpng')

No value below threshold found. Will use first local minimum instead
No value below threshold found. Will use first local minimum instead
xyz: tau = 15.3333

Estimate the embedding dimension (Figure 2b in the article)

[fnnPercent, embeddingDimension] = mdFnn(data, round(tau));
set(gca,'FontSize',fontSize,'fontWeight','normal')
print('Figure2b','-dpng')

Alternative method to find time delay using first local minimum criterion

tau = mdDelay(data, 'maxLag', 25, 'plottype', 'all', 'criterion', 'localMin');
set(gca,'FontSize',fontSize,'fontWeight','normal')
disp('xyz: tau = ' + string(tau))

xyz: tau = 16.6667

Plot the average AMI and standard deviation

tau = mdDelay(data, 'maxLag', 25, 'plottype', 'mean');
set(gca,'FontSize',fontSize,'fontWeight','normal')
disp('xyz: tau = ' + string(tau))

No value below threshold found. Will use first local minimum instead
No value below threshold found. Will use first local minimum instead
xyz: tau = 15.3333

Time delay for the y-variable

tau = mdDelay(data(:,2), 'maxLag', 25, 'plottype', 'all');
set(gca,'FontSize',fontSize,'fontWeight','normal')
disp('y: tau = ' + string(tau))

No value below threshold found. Will use first local minimum instead
y: tau = 15

Time delay for the z-variable

tau = mdDelay(data(:,3), 'maxLag', 25, 'plottype', 'all');
set(gca,'FontSize',fontSize,'fontWeight','normal')
disp('z: tau = ' + string(tau))

z: tau = 12

Time delay for the x and y variables

tau = mdDelay(data(:,1:2), 'maxLag', 25, 'plottype', 'all');
set(gca,'FontSize',fontSize,'fontWeight','normal')
disp('xy: tau = ' + string(tau))

No value below threshold found. Will use first local minimum instead
No value below threshold found. Will use first local minimum instead
xy: tau = 17

Time delay for the x and z variables

tau = mdDelay(data(:,[1,3]), 'maxLag', 25, 'plottype', 'all');
set(gca,'FontSize',fontSize,'fontWeight','normal')
disp('xz: tau = ' + string(tau))

No value below threshold found. Will use first local minimum instead
xz: tau = 15.5

Time delay for the y and z variables

tau = mdDelay(data(:,2:3), 'maxLag', 25, 'plottype', 'all');
set(gca,'FontSize',fontSize,'fontWeight','normal')
disp('yz: tau = ' + string(tau))

No value below threshold found. Will use first local minimum instead
yz: tau = 13.5

Published with MATLAB® R2017b