juliadynamics/transitionsintimeseries.jl

Transition Indicators / Early Warning Signals / Regime Shifts / Change Point Detection

TransitionsInTimeseries.jl

TransitionsInTimeseries.jl is a free and open-source software to easily analyse transitions within timeseries in a reproducible, performant, extensible and reliable way. In contrast to other existing software with similar target application, TransitionsInTimeseries.jl defines a generic interface for how to find transitions and how to test for significance. Within this interface, it is easy to expand the software in three orthogonal ways:

1. Provide the analysis pipelines with new indicators, which can be either self-written or imported from other packages. In particular, the latter offers thousands of metrics that can indicate transitions right out of the box.
2. Add new analysis pipelines for finding transitions.
3. Add new ways for significance testing.

TransitionsInTimeseries is a registered Julia package and can be installed by running:

] add TransitionsInTimeseries

All further information is provided in the documentation, which you can either find online or build locally by running the docs/make.jl file.

Alternative names for this package could have been: Early Warning Signals / Resilience Indicators / Regime-Shift Identifiers / Change-Point Detectors, or however else you want to call them!

transitionsintimeseries.jl's Issues

Collect and list more "Change Point Detection" methods.

Apparently, in the literature of nonlinear dynamics, there are two nomenclatures on regime shifts. One call them "regime transitions" or "regime shifts" or "tipping points". The other calls them "Change Points Detection". It's the same thing. But we were missing a lot of literature because, at least I, was not aware of the second nomenclature.

The following paper is a review on "change points".

A survey of methods for time series change point detection by Samaneh Aminikhanghahi & Diane J. Cook in Knowledge and Information Systems volume 51, pages339–367 (2017)

We should have a look, and little by little add methods from that paper to the list of issues here that request new methods (open new issue -> "new transition indicator").

Neural networks as TIs

Indicator summary

Some pre-trained artificial neural network have shown to be able to predict transitions (and their underlying normal form). I suppose it works well if data corresponds to training set (relatively large for the below-mentioned reference), worse if it is not the case.

Reference

Bury et al. (2021), Deep learning for early warning signals of tipping points.

Codebase

Rely on some deep-learning library in julia (Flux.jl standard but canbe discussed)

Implementation plan

Find and incorporate pre-trained network.

Integrate/incorporate https://github.com/STOR-i/Changepoints.jl

opening this quickly before I forget, and will write more later. Software with similar goals exists and we have to think how to integrate or incorporate to not duplicate effort or content: https://github.com/STOR-i/Changepoints.jl

Register package

To register the package we need to fix the complains in JuliaRegistries/General#93396 and then tag the Julia registrator.

RidgeRegression is clunky to use in new interface

In our new interface we do:

indicators = [var, ar1_whitenoise]
ind_conf = IndicatorsConfig(indicators;
    width = 400
)
# use spearman correlation for both indicators
change_metric = spearman
sig_conf = SignificanceConfig(change_metric;
    width = 20, n_surrogates = 10_000,
    surrogate_method = RandomFourier(),
)

# perform the full analysis
result = indicators_analysis(input, ind_conf, sig_conf)

Using the new RidgeRegression in this new interface is super clunky because we need to deduce both the time vector, and the width of the final change metric timeseries, which the user does not. The user only picks the window. We need to have an internal step where we initialize a RidgeRegression with 0 arguments, and once the function indicators_analysis starts, it optimizes the computation once at the start and then re-uses the type later on.

Use hypothesis testing from TimeseriesSurrogates.jl

Refactor significance testing with: https://juliadynamics.github.io/TimeseriesSurrogates.jl/stable/#Hypothesis-testing

update readme

the package is registered now so at least this part of the readme is false.

Benchmarks in the paper should be done using benchmarking packages

Why are you using the now() - prevtime method for benchmarking and are not using BenchmarkTools and timeit in python?

In the current setup of the julia benchmarks the first compilation time is included in the first run.

openjournals/joss-reviews#6464

New IndicatorChangesConfig: Continuous Linear Segments

This issue is a continuation of the discussion originally in STOR-i/Changepoints.jl#36 .

There we had the quote:

From my point of view, is the very important part of "change analysis" the continuous piece-wise linear regression. See for info the following paper. The method presented in this paper is similar to the standard PELT method for linear segments estimation, but it is moreover forced to be continuous (without step changes between linear segments!!!). See the code on GitLab, too.

The paper is "A new and general approach to signal denoising and eye movement classification based on segmented linear regression". It appears to be a blend between standard PELT and estimating changes in slope in indicator timeseries which @JanJereczek is interested in.

Following the outline I gave in #52 , one could define yet another analysis pipeline type that does the "PELT-like-continuus-linear-segments".

Handle unevenly sampled time series

Is your feature request related to a problem? Please describe.

Empirical data are often unevenly sampled. In the state of the art, an interpolation is often performed but this might bias the result.

In TransitionIndicators.jl we want to offer both options: (1) the user can interpolate their data and pass it or (2) directly pass their unevenly sampled time series, in which case no interpolation happens and all the difficulties are handled internally. (Option 2) might have reduced features compared to (Option 1) but introduces no error.

We won't focus too much on (Option 2) for now but defining the basics for the API to handle it should be a step in near future.

Describe the solution you'd like

Here goes an idea of solution. It does not rely on defining an iterator (maybe a downside?) but has a quite sparse syntax:

using TransitionIndicators

abstract type TimeSeries end

struct EvenspacedTS <: TimeSeries
    nt::Int64
    t::AbstractVector
    spacedim::Int64
    x::AbstractArray
    dt::Real
end

struct UnevenspacedTS <: TimeSeries
    nt::Int64
    t::AbstractVector
    spacedim::Int64
    x::AbstractArray
end

function TimeSeries(t, x)
    nt = length(t)
    spacedim = ndims(x) - 1
    if isequispaced(t)
        dt = mean(diff(t))
        return EvenspacedTS(nt, t, spacedim, x, dt)
    else
        println(
            "\n" *
            " Unevenly spaced time series are work in progress. " * 
            "If you spot bugs, please open an issue on GitHub!\n " * 
            "Bare in mind that unevenly spaced time series come with limited functionalities! \n"
        )
        return UnevenspacedTS(nt, t, spacedim, x)
    end
end

struct WindowView
    width::Int
    stride::Int
    windows::Vector{UnitRange{Int}}
end

function window_view(x, wv::WindowView)
    return [view(x, window) for window in wv.windows]
end

function WindowView(ts::EvenspacedTS; width::Int = default_window_width(ts), stride::Int = default_window_stride(ts))
    si = width:stride:ts.nt
    windows = [k-width+1:k for k in si]
    return WindowView(width, stride, windows)
end
default_window_width(ts::EvenspacedTS) = ts.nt ÷ 100
default_window_stride(ts::EvenspacedTS) = 1

function WindowView(ts::UnevenspacedTS; twidth::Real = default_window_width(ts), stride::Int = default_window_stride(ts))
    ibegin = find_nearest(twidth, ts.t)
    if ibegin < 1
        ibegin = 1
    end
    windows = [findspan(t, twidth, ts.t) for t in ts.t[ibegin:stride:end]]
    return WindowView(0, stride, windows)
end
default_window_width(ts::UnevenspacedTS) = last(ts.t) / 100
default_window_stride(ts::UnevenspacedTS) = 1

function findspan(t, twidth, t_uneven)
    i1 = find_nearest(t-twidth, t_uneven)
    i2 = find_nearest(t, t_uneven)
    if i1 == i2
        error("Time series not dense enough before t=$t for the chosen window width")
    else
        return i1:i2
    end
end

function find_nearest(t, tvec::AbstractVector)
    return argmin((tvec .- t).^2)
end

t = 0:0.1:100
x = sin.(t)
ts = TimeSeries(t, x)
wv = WindowView(ts)
y = zeros(eltype(x), length(wv.windows))
@btime map!(var, y, window_view(x, wv))

tu = cumsum(0.2 .* rand(length(t)))
xu = sin.(t)
tsu = TimeSeries(tu, xu)
wvu = WindowView(tsu)
yu = zeros(eltype(x), length(wvu.windows))
@btime map!(var, yu, window_view(xu, wvu))

Describe alternatives you've considered

Additional context

One open question is whether the window width on unevenly time series should be defined by the number of points or by a time span. I think the latter one is more rigorous.

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Network representations of attractors for change point detection

Indicator summary

TODO; I just got cited by the paper, still gotta read it!

Abstract:

A common approach to monitoring the status of physical and biological systems is through the regular measurement of various system parameters. Changes in a system’s underlying dynamics manifest as changes in the behaviour of the observed time series. For example, the transition from healthy cardiac activity to ventricular fibrillation results in erratic dynamics in measured electrocardiogram (ECG) signals. Identifying these transitions—change point detection—can be valuable in preparing responses to mitigate the effects of undesirable system changes. Here, we present a data-driven method of detecting change points using a phase space approach. Delay embedded trajectories are used to construct an ‘attractor network’, a discrete Markov-chain representation of the system’s attractor. Once constructed, the attractor network is used to assess the level of surprise of future observations where unusual movements in phase space are assigned high surprise scores. Persistent high surprise scores indicate deviations from the attractor and are used to infer change points. Using our approach, we find that the attractor network is effective in automatically detecting the onset of ventricular fibrillation (VF) from observed ECG data. We also test the flexibility of our method on artificial data sets and demonstrate its ability to distinguish between normal and surrogate time series.

Reference

https://www.nature.com/articles/s42005-023-01463-y

Codebase

https://github.com/eugenetkj98/AttractorNetworksPublic

Implementation plan

TODO, haven't read the paper yet.

cc @eugenetkj98 if they are interested to participate here, perhaps they can already provide the implementation strategy.

Eigenvalues of the covariance matrix as early warning signals for critical transitions in ecological systems

Indicator summary

paper abstract:

Many ecological systems are subject critical transitions, which are abrupt changes to contrasting states
triggered by small changes in some key component of the system. Temporal early warning signals such
as the variance of a time series, and spatial early warning signals such as the spatial correlation in a
snapshot of the system’s state, have been proposed to forecast critical transitions. However, temporal
early warning signals do not take the spatial pattern into account, and past spatial indicators only
examine one snapshot at a time. In this study, we propose the use of eigenvalues of the covariance
matrix of multiple time series as early warning signals. We first show theoretically why these indicators
may increase as the system moves closer to the critical transition. Then, we apply the method to
simulated data from several spatial ecological models to demonstrate the method’s applicability.
This method has the advantage that it takes into account only the fluctuations of the system about
its equilibrium, thus eliminating the effects of any change in equilibrium values. The eigenvector
associated with the largest eigenvalue of the covariance matrix is helpful for identifying the regions that
are most vulnerable to the critical transition.

Reference

Scientific Reports | (2019) 9:2572 | https://doi.org/10.1038/s41598-019-38961-5

Shiyang Chen1, Eamon B. O’Dea2,3, John M. Drake2,3 & Bogdan I. Epureanu1

Outline for estimating change points / tipping points

Alright, so we have some measures that have been shown useful to detect transitions. We have some stuff from Early Warning Signals, where an estimator increases more and more towards the thransition, such as AR1-coefficient or Variance. We also have stuff from traditional nonlinear TSA that detect changes in dynamic behavior, such as permutation entropy.

Detecting a change function

As far as I can tell we need two different functions: one that finds significant change in an estimator, and one that finds a continuous increase of an estimator; these are two different things. So we could have: detect_change_points(x, measure, args...; kwargs...) and detect_increase(x, measure, args...; kwargs...).

Estimating significance of results

It's cool that we get the change points, how do you estimate significance? One way is with surrogate tests. E.g., for the slope / early warning signals: you collect the slopes of the surrogates. So, for each rolling window we do a surrogate test. In the surrogate test we collect a "significance" number: this is the normalized distance from the distribution: the absolute distance of the measured value of the slope at this window, minus the mean of the distribution, divided by the width of the central 5-95% quantile. If this number is less than 1 we have no significant result.

(So, this is another function on top of detect_increase or whatever)

sgen = surrogenerator(x, Method)

for i in 1:total_surrogates
   s = sgen()
   ar1_surrogate = sliding_window_apply(ar1_coeff, s, window_viewer)
   slopes = extract_slopes(ar1_surrogates
   matrix_of_slopes[i, :] .= slopes
end

significances = zeros(length(ar1_slopes_real))
for j in 1:size(matrix_of_slopes, 2)
   quantile_595 = quantile(matrix_of_slopes[:, j])
   significances[j] = abs(real_slope - mean(matrix_of_slopes[:, j])/quantile
end

pvalues counting for surrogate significance is broken

Have a look here:

https://juliadynamics.github.io/TransitionsInTimeseries.jl/dev/examples/logistic/#Significance-via-random-Fourier-surrogates

this exacmple used to be perfectly fine. Now, every single time point gets a significance flag. I believe this is because of our re-worked p-value counting code. That code is also over-engineered. We should simplify it back to its original, most simple possible version, because the performance of this code does not matter at all: it has miniscule impact on overall runtime.

New IndicatorsChangesConfig: PELTSegmentation

Alright, with the v0.1 interface out we have a clear way of adding new analysis pipelines for estimating transitions in timeseries. We have the supertype IndicatorsChangesConfig. New subtypes need to extend the estimate_indicator_changes function.

I propose here a new subtype called something like PELTSegmentation that utilizes the existing https://github.com/STOR-i/Changepoints.jl package.

The idea is as follows:

1. The user provides an indicator, and a sliding window configuration. This produces indicator timeseries. Alternatively, the user may give identity as an indicator in which case this first step is skipped.
2. The PELT algorithm from ChangePoints.jl is applied directly to the indicator timeseries. The user provides the parameters of the PELT algorithm as arguments to the new type PELTSegmentation.
3. The results type will then contain the indicator timeseries as well as the time flags that the PELT algorithm identified as change points.

Testing for significance: the approach of surrogates we have currently does not work out of the box, because in the end we will have discrete time points where there are transitions. That's how PELT works: it doesn't give you a change metric (the PELT algorithm integrates the changes of distributions into the resulting optimized change points). However, we can still use surrogates! We can perform the same 1-2-3 step analysis for 1000 surrogates. For each, we get the potential change points. We then have the user specify a significance window: a time window around the change points of the original timeseries. We can count how many surrogates are within this window and if there are many then the transition could happen by chance.

Traceability and dynamical resistance of precursor of extreme events

Indicator summary

The idea is based on developing an indicator that measures "dynamical coupling".

From the abstract:

. Here we introduce a novel time-series- based and non-perturbative approach to efficiently monitor dynamical resistance and apply it to high- resolution observations of brain activities from 43 subjects with uncontrollable epileptic seizures. We gain surprising insights into pre-seizure dynamical resistance of brains that also provide important clues for success or failure of measures for seizure prevention. The novel resistance monitoring perspective advances our understanding of precursor dynamics in complex spatio-temporal systems with potential applications in refining control strategies.

Reference

http://dx.doi.org/10.1038/s41598-018-38372-y

Codebase

The paper does not share code, but I've met Klaus (last author) personally a couple of times and maybe I could ask him to ask his student whether they are willing to share the code.

Optimize codebase

We've spent a lot of time devising the new interface, but practically 0 seconds optimizing it. There are memory "leakages" all over the place. Running:

julia> @btime indicators_analysis($input, $ind_conf, $sig_conf);
  4.146 s (20828885 allocations: 3.07 GiB)

This amount of allocations is of course nonsense. There are probably a lot of incorrect allocations done in windowmap! which shouldnt happen because we use the in-place version with pre-allocated dummies.

Generally speaking, almost all parts of the code base can be optimized. We can start by detecting the bottlenecks. Using the profiling functionality of VSCode we can get a flame graph of what parts of the code are the slowest. Then, we can benchmark the slow functions in separation and see what goes wrong. Do we have type instabilities? Do we allocate unecessary vectors?

Another point to improve is that our functions are not written in a way that allows for easy optimizations. Normally you would want to have a setup function that initializes all containers, and then a function that does the compitations without initializing anything.

Unify metrics computation for original and surrogate timeseries

So far, estimate_indicator_changes and significant_transitions have very similar code but don't call the same, easier to maintain, functions. We should unify this such that both of them call something like:

initialize_indicator_changes()
estimate_indicator_changes!()

Code example in JOSS paper does not work

Trying the code example in the JOSS paper I first got an residual not defined error, but when I replace residual with data I get the following error:

julia> flags = significant_transitions(results, signif)
ERROR: MethodError: no method matching length(::Symbol)

Closest candidates are:
  length(::Combinatorics.Partition)
   @ Combinatorics ~/.julia/packages/Combinatorics/Udg6X/src/youngdiagrams.jl:8
  length(::DataStructures.EnumerateAll)
   @ DataStructures ~/.julia/packages/DataStructures/95DJa/src/multi_dict.jl:96
  length(::Pkg.Types.Manifest)
   @ Pkg ~/.julia/juliaup/julia-1.10.2+0.x64.linux.gnu/share/julia/stdlib/v1.10/Pkg/src/Types.jl:315
  ...

Stacktrace:
 [1] sanitycheck_tail(tail::Symbol, n_ind::Int64)
   @ TransitionsInTimeseries ~/Documents/TransformationsInTimeseriesReview/dev/TransitionsInTimeseries/src/significance/surrogates_significance.jl:177
 [2] significant_transitions(res::SegmentedWindowResults{…}, signif::SurrogatesSignificance{…})
   @ TransitionsInTimeseries ~/Documents/TransformationsInTimeseriesReview/dev/TransitionsInTimeseries/src/significance/surrogates_significance.jl:120
 [3] top-level scope
   @ REPL[141]:1
Some type information was truncated. Use `show(err)` to see complete types.

This is part of the JOSS review openjournals/joss-reviews#6464

definition of equispaced_step for OffsetArray errors

I am going through the untested functions in codecov and the definition of equispaced_step looks wrong for me for OffsetVectors and it throws an error:

using OffsetArrays, TransitionsInTimeseries

julia> a = OffsetVector([-4, 2:10...], -1:8)
julia> equispaced_step(a)
ERROR: BoundsError: attempt to access 10-element OffsetArray(::Vector{Int64}, -1:8) with eltype Int64 with indices -1:8 at index [9]
Stacktrace:
 [1] throw_boundserror(A::OffsetVector{Int64, Vector{Int64}}, I::Tuple{Int64})
   @ Base ./abstractarray.jl:737
 [2] checkbounds
   @ ./abstractarray.jl:702 [inlined]
 [3] getindex
   @ ~/.julia/packages/OffsetArrays/E3eDu/src/OffsetArrays.jl:437 [inlined]
 [4] equispaced_step(t::OffsetVector{Int64, Vector{Int64}})
   @ TransitionsInTimeseries ~/.julia/packages/TransitionsInTimeseries/XwpWT/src/misc/timeseries.jl:24
 [5] top-level scope
   @ REPL[9]:1

Also this function should be tested.

This is part of the JOSS review openjournals/joss-reviews#6464

What's whitenoise in AR1?

I don't understand why there is _whitenoise at the end of the function ar1_whitenoise. The docstirng also does not describe it. We need to either explain in the docstirng or simplify the function to just ar1.

Improve tests: better tests and more coverage

At the moment, the tests of this package are in an poor state. I think that is normal. The package is in early stages of development and hence developer time is spent mainly building the library API and building a sufficient pool of features. Nevertheless, as the API stabilizes, we should shift development effort more into improving tests.

First, we need to write much more tests and improve our test coverage (i.e., the total amount of functions from the source code that are called and tested in the test suite). Second, we should write good quality tests. Here is a guideline for good tests:

• Actually unit: test atomic, self-contained functions. Each test must test only one thing, the unit. When a test fails, it should pinpoint the location of the problem. Testing entire processing pipelines (a.k.a. integration tests) should be done only after units are covered, and only if resources/time allow for it!
• Known output / Deterministic: tests defined through minimal examples that their result is known analytically are the best tests you can have! If random number generation is necessary, either test valid output range, or use seed for RNG
• Robust: Test that the expected outcome is met, not the implementation details. Test that the target functionality is met without utilizing knowledge about the internals. Also, never use internal functions in the test suite.
• High coverage: the more functionality of the code is tested, the better
• Clean slate: each test file should be runnable by itself, and not rely on previous test files
• Fast: use the minimal amount of computations to test what is necessary
• Regression: Whenever a bug is fixed, a test is added for this case
• Input variety: attempt to cover a wide gambit of input types

One doesn't have to worry about re-writing all tests. In fact, a PR "correcting" a single test file is already very much welcomed!

Improve the performance of `PrecomputedLowfreqPowerSpectrum`

Indicator summary

Compute ratio between low-frequency power spectrum and total power spectrum. What is "low-frequency"? Defined by user but standard is something like 10 lowest percent of the total power spectrum. Performs well wherever CSD is expected (conversely bad wherever CSD might be "hidden", e.g. by high signal-to-noise ratio).

Reference

Bury et al. (2020), Detecting and distinguishing tipping points using spectral early warning signals.

Codebase

Rely on FFTW for FFT

Implementation plan

struct LowfreqPowerSpectrum <: Function
    plan::Any
    i_lofreq::Int
end

function LowfreqPowerSpectrum(x::AbstractVector, width::Int; q_lofreq = 0.1)
    p = plan_rfft(x[1:width+1])
    fft_x = p*x
    i_lofreq = Int(round(length(fft_x) * q_lofreq))
    return LowfreqPowerSpectrum(p, i_lofreq)
end

function (lfps::LowfreqPowerSpectrum)(x::AbstractVector)
    ps = abs.(lfps.plan * x)
    return sum(ps[1:lfps.i_lofreq])
end

AR1 for non-stationary correlated noise

Indicator summary

Correlated noise can lead to bad transition recognition if standard CSD is applied (variance, lag1-autocorelation, AR1 parameter). Therefore use an adapted formula that addresses this. Works bad for small window sizes.

Reference

• Boettner and Boers (2021), Critical Slowing Down In Dynamical Systems Driven By Non-Stationary Correlated Noise

Also of interest in that context: work of Frank Kwasniok (in preparation?)

Codebase

Nope

Implementation plan

See eq (17-20) of Boettner and Boers (2021)

Write example for "Automatic detection of abrupt transitions in paleoclimate records"

Paper: https://aip.scitation.org/doi/10.1063/5.0062543

Now that the basic interface for TransitionsIndicators.jl is finished, we can re-create this paper. What we have to do is use identity as indicator and write a trivial function that re-uses the existing nonparametric Kolmogorov–Smirnov from HypothesisTests.jl as the change metric.

Kolmogorov-Smirnov test as an indicator of a transition

Indicator summary

From the abstract;

Here, we present a robust methodology for detecting abrupt transitions in proxy records that is applied to ice core and speleothem records of the last climate cycle. This methodology is based on the nonparametric Kolmogorov–Smirnov (KS) test for the equality, or not, of the probability distributions associated with two samples drawn from a time series, before and after any potential jump. To improve the detection of abrupt transitions in proxy records, the KS test is augmented by several other criteria and it is compared with recurrence analysis. The augmented KS test results show substantial skill when compared with more subjective criteria for jump detection. This test can also usefully complement recurrence analysis and improve upon certain aspects of its results.

This method has been used to generate the PaleoJump library: https://paleojump.github.io/

Reference

https://aip.scitation.org/doi/abs/10.1063/5.0062543

https://arxiv.org/abs/2206.06832

Codebase

Nope, but maybe @paleojump is willing to share the code base.

Implementation plan

Once we have the Kolmogorov-Smirnov test implemented, the rest should be straight forward given that the moving window functionality is already implemented by Jan.

Delayed slowing down

Indicator summary

Matthias Marconi et al has shown that critical slowing down are not good indicators for systems with a parameter changing in time. They made a new one based on their findings which they call "delayed slowed down".

I came across this work while attending a presentation in Aberdeen Dynamics Days 2022.

Reference

https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.125.134102

Arxiv version exists as well. Although, I am not sure whetheer the paper has all the recent results that the presentation showed. We can ask the lead author to share newer paper (will do so).

Codebase

No.

Implementation plan

Not yet.

Increase test coverage to at least 85%

Local /finite-time Lyapunov exponents as early warning signals

Indicator summary

Abrupt climate change has an important impact on sustainable economic and social development, as well as ecosystem.
However, it is very difficult to predict abrupt climate changes because the climate system is a complex and nonlinear system.
In the present paper, the nonlinear local Lyapunov exponent (NLLE) is proposed as a new early warning signal for
an abrupt climate change. The performance of NLLE as an early warning signal is first verified by those simulated abrupt
changes based on four folding models. That is, NLLE in all experiments showed an almost monotonous increasing trend as
a dynamic system approached its tipping point. For a well-studied abrupt climate change in North Pacific in 1976/1977, it
is also found that NLLE shows an almost monotonous increasing trend since 1970 which give up to 6 years warning before
the abrupt climate change. The limit of the predictability for a nonlinear dynamic system can be quantitatively estimated by
NLLE, and lager NLLE of the system means less predictability. Therefore, the decreasing predictability may be an effective
precursor indicator for abrupt climate change

Reference

Climate Dynamics (2021) 56:3899–3908
https://doi.org/10.1007/s00382-021-05676-1
Decreasing predictability as a precursor indicator for abrupt climate
change
Wenping He et al

Implementation plan

To be discussed. Local/nonlinear/finite-time Lyapunov Exponents are implemented in ChaosTools.jl but for analytic dynamical systems, not timeseries....

Collection of some Early Warning Signals / Resilience Indicators / Regime Shift Identifiers / Transition Identifiers

EDIT: These algorithms will live in https://github.com/JuliaDynamics/TransitionIdentifiers.jl

---

In the near future we'll be adding a module in DynamicalSystems.jl about functionality related to "early warning signals, resilience indicators, tipping indicators, precurson signals, quantities that identify regime shifts" and similar functionality that is typically used in tipping points analysis. The plan is to have both timeseries analysis tools, but also dynamical-system-based tools (similarly to how you can compute a lyapunov exponent from a dynamical system / ODE with lyapunov or use a reconstructed signal and get the exponent from data with lyapunov_from_data).

The purpose of this Issue is to collect existing functionality online and in the current literature, and serve both as a starting point, but also as an indicator of what new research directions are useful to be pursued in the future.

(will be editing frequently this issue to add more stuff)

Methodologies/Methods/Quantities

• Physical precursor signals. Such signals are context-specific, and depend on the physics and causal connections of the signal at hand. They can be utilized to anticipate (TODO: write more).
• Critical Slowing Down indicators. In the specific case of a sadle-node bifurcation in a 1D (or, quasi-1D-potential) case driven by weak, Gaussian noise, it has been showing that the autocorrelation and variance of a signal (computed via a sliding window) increase the closer the system gets to the bifurcation. TODO; Add references.
• Kolmogorov-Smirnov difference. TODO: write more, Witold's work.
• Recurrence quantificaiton analysis. TODO: write more (also Witold's work).
• basins fractions. TODO: Write more (Menck2013 and our work with Alex).

Existing code bases

• Vasilis Dakos toolbox/software TODO write more

Covariant Lyapunov vectors as indicators

Indicator summary

I came across this from a presentation of Sarah Hallerberg in Dynamics Days 2022 Aberdeen. They worked on how covariant lyapunov vectors behave close to a critical transition. I quote from their abstract:

Studying critical transitions in several models of fast-slow systems, i.e., a network of coupled FitzHugh-Nagumo oscillators, models for Josephson junctions, and the Hindmarsh-Rose model, we find that tangencies between covariant Lyapunov vectors are a common and maybe generic feature during critical transitions.

more importantly:

In the presence of noise, we find the alignment of covariant Lyapunov vectors and changes in finite-time Lyapunov exponents to be more successful in announcing critical transitions than common indicator variables as, e.g., finite-time estimates of the variance.

Reference

https://journals.aps.org/pre/abstract/10.1103/PhysRevE.96.032220

and also this paper for the method to estimate CLVs from data: https://paperswithcode.com/paper/estimating-covariant-lyapunov-vectors-from

Codebase

Nope. But I can ask Nahal Sharafi, perhaps she is willing to share.

Implementation plan

Not yet.

Automatic flag creation function

Currently we do:

signif_idxs = vec(count(results.pval .< threshold, dims = 2) .>= 2)

to find the flags of transitions. This should be automated to affunction transition_flags that takes as an argument the results type and a second argument for which indicators to detect flags

least square density difference for "change point" detection

Indicator summary

Alright, I just saw this Julia code, https://github.com/johncwok/ChangePointDetection.jl , that implements an algorithm called least square density difference (LSDD) method. It is used to detect changepoints in time-series or to infer wether or not two time-series come from the same underlying probability distribution.

Reference

The LSDD method was developped in the article Density-difference estimation from M. Sugiyama, T. Suzuki, T. Kanamori, M. C. du Plessis, S. Liu, and I. Takeuchi. in 2013.

http://www.ms.k.u-tokyo.ac.jp/sugi/2013/LSDD.pdf

Codebase

https://github.com/johncwok/ChangePointDetection.jl

Implementation plan

Best to have a look at the code and see if it can be improved performance wise. Probably also integrated with our sliding window viewer.

Detrended Fluctuation Analysis as an early warning indicator

Indicator summary

We assess the proximity of a system to a bifurcation point using a degenerate fingerprinting method that estimates the declining decay rate of fluctuations in a time series as an indicator of approaching a critical state. The method is modified by employing Detrended Fluctuation Analysis (DFA) which improves the estimation of short-term decay, especially in climate records which generally possess power-law correlations. When the modified method is applied to GENIE-1 model output that simulates collapse of the Atlantic thermohaline circulation, the bifurcation point is correctly anticipated. In Greenland ice core paleotemperature data, for which the conventional degenerate fingerprinting is not applicable due to the short length of the series, the modified method detects the transition from glacial to interglacial conditions. The technique could in principle be used to anticipate future bifurcations in the climate system, but this will require high-resolution time series of the relevant data.

Reference

A modified method for detecting incipient bifurcations in a dynamical system
V. N. Livina, T. M. Lenton

GEOPHYSICAL RESEARCH LETTERS, VOL. 34, L03712, doi:10.1029/2006GL028672, 2007

https://doi.org/10.1029/2006GL028672

Implementation plan

There is DFA code in Julia, but the repo doesn't look super trustworthy: https://github.com/johncwok/DetrendedFluctuationAnalysis.jl (hasn't been updated in 4 years). I know the author, and I can ask if the repo is LGTM still. I've used used DFA extensively in the past in the context of music timeseries analysis, I'm sure this one will be easy to tackle.

Scale-averaged wavelet coefficient and local Hurst exponent

Indicator summary

Scale-averaged wavelet coeff is used to estimate the variance of high-frequency fluctuations and local Hurst exponent is interpreted as an estimator of the correlation time in the high frequencies. Here in same issue because published in same paper and might be specifically suited for DO-events. Both are CSD indicators and hence come with usual caution.

Reference

Boers (2018), Early-warning signals for Dansgaard-Oeschger events in a high-resolution ice core record

Codebase

Not that I would know

Implementation plan

Nothing so far

Cannot reproduce benchmarks ewstools_tuto-1.csv missing

I can not reproduce the benchmarks of TransitionsInTimeseries.jl because the file ewstools_tuto-1.csv is missing.

Could you please provide this file in the paper/code folder?
This is part of the JOSS review openjournals/joss-reviews#6464

Restoring rate

Indicator summary

Linearizing a 1D nonlinear model, one obtains $\dot{x} = \lambda , x$, with $\lambda$ the restoring rate. It can be estimated in a way that is robust to changes in variance of the noise process leading to the fluctuations around the attractor. The restoring rate is a CSD indicator

Reference

Boer et al. (2021), Observation-based early-warning signals for a collapse of the Atlantic Meridional Overturning Circulation.

Codebase

Not that I would know but should be easy to implement.

Implementation plan

C.f. Boers et al. (2021), "Methods" section

Convenience functions for plotting

Is your feature request related to a problem? Please describe.
Plotting of the results is generic. Therefore providing a convenience function for it would cut down the amount of code needed.

Describe the solution you'd like
Rely on CairoMakie and allow something like:

plot(results, signif)

Add developers docs page: Adding new "anything"

Due to our design, adding new "something" has well defined steps: new indicator, change metric, estimation of changes, estimation of significance.

We should have a "developers docs" pages that outlines what are the steps necessary to add a new method. This will further enhance the claims we make in the paper about accessiblity.