Classical RQA measures for periodic and chaotic orbits

[Datseris]

Here are two parameters for the Lorenz system:

One is clearly periodic, the other is clearly chaotic. I now want to demonstrate the use of the classical RQA measures in these two cases:

lor = Systems.lorenz()
for (i, ρ) in enumerate((69.75, 28.0))
    println("\nρ = $ρ ")
    set_parameter!(lor, 2, ρ)
    ε = 2.0
    t, dt = 20.0, 0.005
    tr = trajectory(lor, t; dt = dt, Ttr = 200.0)
    R = RecurrenceMatrix(tr, ε)

    println(determinism(R))
    println(dl_max(R))
    println(dl_average(R))
    println(dl_entropy(R))
end

ρ = 69.75 
0.9996688958227877
2598
83.35871404399323
4.626554662807655

ρ = 28.0 
0.9997196293101414
4000
74.48904006046862
4.938120518819617

Here are my questions:

1. Almost all parameters are identical between the measures. How can this be / How can I understand this? I expected big differences at least in the entropies.
2. The maximum diagonal for thechaotic case is higher than the regular case. Conceptually I cannot imagine this being true. Why does it happen?
3. Can I actually use the classical RQA measures to point out differences between periodic and (simple) chaotic trajectories, or should I go about a different approach in the video tutorial?

[heliosdrm]

There are various things here, related to parameter selection:

1. The thresholds used in the "periodic" and "chaotic" versions of the Lorenz system are not equivalent. The phase space of the periodic system is about two times the other (see the images below, and the values of distancematrix(tr)), so it would make more sense to use a threshold of 2ε in the periodic version. Actually, the fact that the longest line in the "periodic" version is around 2,500 means with the current thresholds you are missing some recurrences in the long diagonal lines (it seems that the longest after the main diagonal is around 3,300). Using 2ε for the periodic version will give you the "correct" (or at least expected) longest diagonal line.

2. The longest line in the chaotic version is 4,000 (the full extent of the trajectory), which is also a warning of a "thick" line around the LOI. I would avise to use theiler=2 (you may use it in both examples). With this and the previous recommendation about the threshold, you will have Lmax=3300 for the periodic, and Lmax=1317 for the chaotic, which look much more like what the plots are telling (let aside the LOI).

The differences in ENTR will also be greater, although not extraordinary. This is because the plot of the periodic version is composed not only of the long diagonal lines corresponding to the main periodicity, but also of shorter tangential fragments of the trajectories, when the "loops" of the system cross to each other -although they were not visible in your plot-:

So, the distribution of diagonal lines is not in the end that much different in both systems. There is a bigger difference in TREND (one or two orders of magnitude, with smaller "fading" in the periodic sytem).

[pucicu]

(1) dl_max: increase the Theiler window in order to avoid counting the very long lines in the direct vicinity of the main diagonal (2) DET is comparable, which makes sense (3) in the RP it is not clearly visible whether the diagonal lines for the periodic case are really continuous or interrupted – they should be continuous and not interrupted; the lines also appear quite think, perhaps reduce ε in order to get thinner lines (preferably of width 1) (4) there is a general problem with the diagonal lines of periodic systems, because they cross the borders of the RP. this results in a “strange” histogram of diagonal line lengths, which is discrete but has entries at different lengths and not, as we would expect, only at one length. this biased histogram leads to biased values of dl_average and dl_entropy, e.g. higher entropy for periodic dynamics than for chaotic. in the strict sense of describing the distribution of lines in the RP this would be correct, but has nothing to do with the complexity of the dynamics of the system.

for issue (4) I’m just thinking whether we should already include now the correction schema which we are currently preparing for a publication. the idea is that all diagonal lines that have their first and last point at the border of the RP (as all the lines of a RP from a periodic system) are removed from the histogram (important: start and end point of the line are at the border, meaning the line is cut twice!). but then these histograms are empty for periodic systems. therefore, the first (longest) line that crosses the borders (and was actually excluded) should be included again. this works fine but this is still work in progress and not published published. some further tests are necessary.

Am 25.01.2019 um 10:29 schrieb George Datseris notifications@github.com:

Here are two parameters for the Lorenz system:

https://user-images.githubusercontent.com/19669089/51736657-ecaf5500-208a-11e9-803e-d496d4a83858.png One is clearly periodic, the other is clearly chaotic. I now want to demonstrate the use of the classical RQA measures in these two cases:

lor = Systems.lorenz() for (i, ρ) in enumerate((69.75, 28.0)) println("\nρ = $ρ ") set_parameter!(lor, 2, ρ) ε = 2.0 t, dt = 20.0, 0.005 tr = trajectory(lor, t; dt = dt, Ttr = 200.0) R = RecurrenceMatrix(tr, ε)

println(determinism(R))
println(dl_max(R))
println(dl_average(R))
println(dl_entropy(R))

end

ρ = 69.75 0.9996688958227877 2598 83.35871404399323 4.626554662807655

ρ = 28.0 0.9997196293101414 4000 74.48904006046862 4.938120518819617 Here are my questions:

Almost all parameters are identical between the measures. How can this be / How can I understand this? I expected big differences at least in the entropies. The maximum diagonal for thechaotic case is higher than the regular case. Conceptually I cannot imagine this being true. Why does it happen? Can I actually use the classical RQA measures to point out differences between periodic and (simple) chaotic trajectories, or should I go about a different approach in the video tutorial? — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/JuliaDynamics/RecurrenceAnalysis.jl/issues/64, or mute the thread https://github.com/notifications/unsubscribe-auth/AA7WkMkaf71u8UbuSFsQ5_d1fC9pCVvyks5vGs6FgaJpZM4aSgZP.

[Datseris]

Thank you both for your comments, I will read them carefully and modify both my understanding and the tutorial!

I only have a quick comment for the last point that Norbert mentioned:

I actually wanted to open a feature request issue for a function that identifies "infinite lines" i.e. lines that both start and end at the borders of the RP. This function should have keyword arguments theiler, border, i.e. same keywords as trend.

This function would allow one to not only see if there is periodic behavior in the system, but also count the period, because the vertical distance between the infinite lines is actually the period of the system (at least an approximation with +/- the sampling rate).

[pucicu]

would be interesting indeed.

Am 25.01.2019 um 14:46 schrieb George Datseris notifications@github.com:

Thank you both for your comments, I will read them carefully and modify both my understanding and the tutorial!

I only have a quick comment for the last point that Norbert mentioned:

I actually wanted to open a feature request issue that identifies "infinite lines" i.e. lines that both start and end at the borders of the RP. This function should have keyword arguments theiler, border, i.e. same keywords as trend.

This function would allow one to not only see if there is periodic behavior in the system, but also count the period, because the vertical distance between then infinite lines is actually the period of the system (at least an approximation with +/- the sampling rate).

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/JuliaDynamics/RecurrenceAnalysis.jl/issues/64#issuecomment-457577302, or mute the thread https://github.com/notifications/unsubscribe-auth/AA7WkEDZjH5b-6_XNdfHHF8WttZsqF-5ks5vGwqkgaJpZM4aSgZP.

[Datseris]

Hi, I am back again with one (hopefully final) question.

For reference, using the parameters you suggested the TREND has the values:

ρ = 69.75
trend: 1.4490199175443764e-6

ρ = 28.0 
trend: 7.85248479721087e-9

it is identical for both, equal to zero.

---

Now, I wanted to make an example that shows that the TREND value clearly increases for a system that is non-stationary. Since the review papers state that TREND will be significantly off zero. So I created a lorenz system where the ρ parameter changes with time. Here is the results:

First, create the system:

function lorenz(u0=[0.0, 10.0, 0.0]; σ = 10.0, ρ = 28.0, β = 8/3)
    return CDS(loop, u0, [σ, ρ, β], loop_jac)
end
function lor_ns(u, p, t)
    @inbounds begin
        σ = p[1]; ρ = u[4]; β = p[2]
        du1 = σ*(u[2]-u[1])
        du2 = u[1]*(ρ-u[3]) - u[2]
        du3 = u[1]*u[2] - β*u[3]
        duρ = 1.0
        return SVector{4}(du1, du2, du3, duρ)
    end
end

nonstationary_lorenz = ContinuousDynamicalSystem(lor_ns, [0, 10, 0, 28.0], [10.0, 8/3])

then, plot a recurrence matrix:

t = 50.0; dt = 0.01
tr = trajectory(nonstationary_lorenz, 50.0; dt = 0.01)
tvec = 0:dt:t
ε = 4.0
R = RecurrenceMatrix(tr[:, 1:3], ε) # notice I use 1:3
x, y = coordinates(R)
figure(figsize = (6,6))
scatter(tvec[x], tvec[y], s = 1, alpha = 0.1, color = "k")
xlim(1, t); ylim(1, t); gca().set_aspect("equal")
xlabel("t"); ylabel("t"); 

Finally compute some RQA parameters:

RQA = rqa(R, theiler = 4, lmin = 2)
println("maximum diagonal: ", RQA.Lmax)
println("average diagonal: ", RQA.L)
println("determinism: ", RQA.DET)
println("trend: ", RQA.TREND)

maximum diagonal: 1044
average diagonal: 25.040637450199203
determinism: 0.9986018430251032
trend: -1.3878710653763777e-5

You can see that besides the fact that L, Lmax and DET are almost identical with previously, the TREND did not change at all, it is still zero.

---

How can I interpret this fact? I was really expecting this not to be the case. I have to say that I have never worked with recurrence plots before starting to make these tutorials, so I am not sure what I am doing wrong. In addition, I am following the reviews in the most naive way possible. I also have no sense of scale for the relevance of the TREND values. Something that is 1e-5 is "zero" for me.

[heliosdrm]

In my experience such small values of TREND are expected. Some authors express TREND values multiplied x10^3 to obtain greater values, e.g. Riley et al. (1999), or Webber & Zbilut (2005) - in both cases the value is actually multiplied x10^5, because they reported percentages of recurrence, not ratios.

I'll do a last PR with such "correction" (without percentage, only x10^3), before the release of 1.0.0. Just merge if you think it is advisable.

[Datseris]

Okay but this still leaves the difference of the nonstationary result to be 10 times larger than the regular/periodic result (and still hover around 0). I was just wondering whether this is a significant difference or not.

I am also wondering whether I should just completely change the example I use for demonstrating recurrence plots; If you have any suggestions please do let me know.

[heliosdrm]

I cannot give much advice about the trend, because it is not reported in as many places as other parameters like RR, DET, ENTR, etc., and I have not used it myself that much. But perhaps this is a significant fact: one of the reasons why I don't use it is that at least in my field it does not seem to be very robust.

For instance, when the signal is not stationary, the length of the measurement affects the TREND a lot in a non-expected way: one might expect that the longer is a nonstationary measurement, the more marked would be the trend, because the signal drifts away and the recurrence plot is nearly empty except in the vicinity of the main diagonal. But the output is the opposite, and longer signals actually have ~higher~ smaller trends!!

The explanation is that TREND is the slope of a linear regression, and in such cases the "paling" effect seen in the plot is not linear at all: recurrences are dense near the main diagonal, and then null in a big area of the plot, such that the bigger is that empty area, the flatter is the slope.

In my opinion a robust TREND parameter should rule out those big empty areas at the ends of the plots (perpendicular to the diagonals) to avoid such false "flattening", but I have not seen any publication commenting on that. - Nevertheless, it may be an interesting warning to present in the tutorial. Your example of the nonstationary Lorenz system is a good one - you only have to change the length of the trajectory to show it!

[Datseris]

Thanks Helios. I've decided to not show TREND at all in my tutorial. From our conversation here I concluded that It simply is not reliable enough. Even in the most basic and straightforward cases, like the examples I have mentioned here, it still can't clearly distinguish between the two obviously different cases of stationarity.

I have changed the timeseries lengths and, as you suggested, indeed it worked; the TREND was much higher for the non-stationary case. I then decided to compare them side by side, but I said "wait a moment, If I am to compare them then I need the same timevector for both". Well, when I did that the process failed, i.e. nonstationary and stationary timeseries TREND in the same order of magnitude (but one was negative and one was positive).

Thus, I don't think TREND is something worth mentioning in a video tutorial that many beginners will be watching. If I say a sentence like "it is an estimate of the stationarity of a timeseries", I don't have any evidence to prove it.

Here is my code, FYI:

function lor_ns(u, p, t)
    @inbounds begin
        σ = p[1]; ρ = u[4]; β = p[2]
        du1 = σ*(u[2]-u[1])
        du2 = u[1]*(ρ-u[3]) - u[2]
        du3 = u[1]*u[2] - β*u[3]
        duρ = 1.0
        return SVector{4}(du1, du2, du3, duρ)
    end
end

nonstationary_lorenz = ContinuousDynamicalSystem(lor_ns, [0, 10, 0, 18.0], [10.0, 8/3])

t = 50.0; dt = 0.05
tr_ns = trajectory(nonstationary_lorenz, t; dt = dt, Ttr = 2.0)

tr_st = trajectory(Systems.lorenz(), t; dt = dt, Ttr = 200.0)
tvec = 0:dt:t
figure(figsize = (10,6))

for (i, tr) in enumerate((tr_st, tr_ns))

    c = i == 1 ? "C3" : "C0"

    x = tr[:, 3]
    subplot2grid((3,2), (0, i-1))
    plot((0:dt:t), x, color = c, lw = 1.0)
    title(i == 1 ? "stationary" : "non-stationary")    

    subplot2grid((3,2), (1, i-1), rowspan = 2)

    ε = 4.0
    R = RecurrenceMatrix(tr[:, 1:3], ε) # notice, I use 1:3 !
    x, y = coordinates(R)
    scatter(tvec[x], tvec[y], s = 1, alpha = 0.1, color = c)
    xlim(0, t); ylim(0, t); gca().set_aspect("equal")
    xlabel("t"); ylabel("t"); 
end

for (i, tr) in enumerate((tr_st, tr_ns))
    println(i == 1 ? "stationary" : "non-stationary")    

    ε = 4.0
    i == 1 && (ε *= 2)
    R = RecurrenceMatrix(tr, ε)

    RQA = rqa(R, theiler = 4, lmin = 4)

    println("maximum diagonal: ", RQA.Lmax)
    println("average diagonal: ", RQA.L)
    println("determinism: ", RQA.DET)
    println("trend: ", RQA.TREND)
    println()
end

stationary
maximum diagonal: 121
average diagonal: 9.019710144927537
determinism: 0.896824024439449
trend: 0.007400889616943176

non-stationary
maximum diagonal: 229
average diagonal: 11.237547892720306
determinism: 0.9049676025917927
trend: -0.01944905911174211

---

I am now thinking of some other intuitive and most importantly robust RQA measure that I could show. I am thinking of using laminarity, as its definition is a simple thing: how much portion of recurrent points form vertical lines.

I was thinking of using the logistic map just before the period 3 window, because I know it is intermittent there. But if you have any other suggestions, of e.g. continuous systems, please do let me know.

[heliosdrm]

In Norbert's paper (https://doi.org/10.1016/j.physrep.2006.11.001) there is a nice example of a "stretched" Rössler system (a=0.15, b=0.2 and c=10, with a transformation of t to t²) that shows curved lines (see snapshot below), where laminarity and other vertical features may be interesting to show.

[Datseris]

Thanks for the tip Helios! I'm using the logistic map now, it turned out okay! With this, I finished the tutorial on the RQA part. You'll be seeing the results soon!

[pucicu]

I agree to not show TREND in your tutorial. Although I repeat myself: in the pitfalls paper in IJBC 2011 you can find some discussion on this measure and how tricky it is. I’m therefore not using it, and obviously, because of the pitfalls it is not really popular.

Am 27.01.2019 um 11:33 schrieb George Datseris notifications@github.com:

Thanks Helios. I've decided to now show TREND at all in my tutorial. From our conversation here I concluded that It simply is not reliable enough. Even in the most basic and straightforward cases, like the examples I have mentioned here, it still can't clearly distinguish between the two obviously different cases of stationarity.

I have changed the timeseries lengths and, as you suggested, indeed it worked; the TREND was much higher for the non-stationary case. I then decided to compare them side by side, but I said "wait a moment, If I am to compare them then I need the same timevector for both". Well, when I did that the process failed, i.e. nonstationary and stationary timeseries TREND in the same order of magnitude (but one was negative and one was positive).

Thus, I don't think TREND is something worth mentioning in a video tutorial that many beginners will be watching. If I say a sentence like "it is an estimate of the stationarity of a timeseries", I don't have any evidence to prove it.

Here is my code, FYI:

function lor_ns(u, p, t) @inbounds begin σ = p[1]; ρ = u[4]; β = p[2] du1 = σ(u[2]-u[1]) du2 = u[1](ρ-u[3]) - u[2] du3 = u[1]u[2] - βu[3] duρ = 1.0 return SVector{4}(du1, du2, du3, duρ) end end

nonstationary_lorenz = ContinuousDynamicalSystem(lor_ns, [0, 10, 0, 18.0], [10.0, 8/3])

t = 50.0; dt = 0.05 tr_ns = trajectory(nonstationary_lorenz, t; dt = dt, Ttr = 2.0)

tr_st = trajectory(Systems.lorenz(), t; dt = dt, Ttr = 200.0) tvec = 0:dt:t figure(figsize = (10,6))

for (i, tr) in enumerate((tr_st, tr_ns))

c = i == 1 ? "C3" : "C0"

x = tr[:, 3]
subplot2grid((3,2), (0, i-1))
plot((0:dt:t), x, color = c, lw = 1.0)
title(i == 1 ? "stationary" : "non-stationary")

subplot2grid((3,2), (1, i-1), rowspan = 2)

ε = 4.0
R = RecurrenceMatrix(tr[:, 1:3], ε) # notice, I use 1:3 !
x, y = coordinates(R)
scatter(tvec[x], tvec[y], s = 1, alpha = 0.1, color = c)
xlim(0, t); ylim(0, t); gca().set_aspect("equal")
xlabel("t"); ylabel("t");

end https://user-images.githubusercontent.com/19669089/51799713-ef45b200-2224-11e9-81ee-8e635ac6cc74.png for (i, tr) in enumerate((tr_st, tr_ns)) println(i == 1 ? "stationary" : "non-stationary")

ε = 4.0
i == 1 && (ε *= 2)
R = RecurrenceMatrix(tr, ε)

RQA = rqa(R, theiler = 4, lmin = 4)

println("maximum diagonal: ", RQA.Lmax)
println("average diagonal: ", RQA.L)
println("determinism: ", RQA.DET)
println("trend: ", RQA.TREND)
println()

end stationary maximum diagonal: 121 average diagonal: 9.019710144927537 determinism: 0.896824024439449 trend: 0.007400889616943176

non-stationary maximum diagonal: 229 average diagonal: 11.237547892720306 determinism: 0.9049676025917927 trend: -0.01944905911174211 I am now thinking of some other intuitive and most importantly robust RQA measure that I could show. I am thinking of using laminarity, as its definition is a simple thing: how much portion of recurrent points form vertical lines.

I was thinking of using the logistic map just before the period 3 window, because I know it is intermittent there. But if you have any other suggestions, of e.g. continuous systems, please do let me know.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/JuliaDynamics/RecurrenceAnalysis.jl/issues/64#issuecomment-457906300, or mute the thread https://github.com/notifications/unsubscribe-auth/AA7WkLhvnznyxmPCF0hfyBgeaLmyvO-Wks5vHYCXgaJpZM4aSgZP.