

NormalisedPredictiveAsymmetryTest

# CausalityTools.CausalityTests.NormalisedPredictiveAsymmetryTest — Type.

NormalisedPredictiveAsymmetryTest(predictive_test::CausalityTest; f::Number = 1.0)

The parameters for a normalised predictive asymmetry causality test [1]. For the non-normalised version, see PredictiveAsymmetryTest.

Mandatory keywords

• predictive_test: An instance of a predictive causality test that explicitly uses prediction lags (e.g. VisitationFrequencyTest or TransferOperatorGridTest).
• f::Number: The normalisation factor.

About the prediction lags

The prediction lags in the predictive causality test must consist of n negative integers and n positive integers that are symmetric around zero.

In other words, negative lags must exactly match the positive lags but with opposite sign. The zero lag can be included, but will be ignored, so it is possible to give ranges too.

Examples

Let's define a a transfer entropy test with which to compute transfer entropies symmetrically around lag zero. We'll use that as input to the predictive asymmetry test.

bin = RectangularBinning(4) # divide each axis into 4 equal-length intervals
ηs = [-5:-1; 1:5] # exclude the zero lag (it is irrelevant for the asymmetry)
test_visitfreq = VisitationFrequencyTest(binning = bin, ηs = ηs)

# Normalising predictive asymmetry values to f = 1.0 times the mean transfer entropy.
NormalisedPredictiveAsymmetryTest(predictive_test = test_visitfreq, f = 1.0)

References

1. Diego, David, Kristian Agasøster Haaga, Jo Brendryen, and Bjarte Hannisdal. A simple test for causal asymmetry in complex systems. In prep.

source

Example

Preparations

First we load necessary packages.

using CausalityTools, DynamicalSystems, Plots

# Use the pyplot backend, so we can use LaTeX formatted strings. To do this,
# you need the PyPlot package installed
using LaTeXStrings; pyplot()

Plots.PyPlotBackend()

For this example, we'll use the built-in logistic4 system, which consists of four unidirectionally coupled logistic maps coupled $x \to y \to z \to w$. With default values, the system blows up for some initial conditions, so we'll start by just making a thin wrapper that finds us a good orbit.

"""
    good_orbit(f::Function, npts::Int; Ttr::Int = 100, max_attempts = 1000)

Draw a good orbit from a dynamical system defined by `f` when called. Useful
when some initial conditions yield attractors while others do not. CAREFUL:
this function is defined recursively.
"""
function good_orbit(f::Function, npts::Int; Ttr::Int = 100, max_attempts = 1000)
    sys = f()
    ts = columns(trajectory(sys, npts, Ttr = Ttr))

    tries = 0
    while tries < max_attempts
        tries += 1
        if all(isfinite.(hcat(ts...)))
            return ts
        end
    end
    println("Did not manage to find good orbits within $(max_attempts) attempts. Re-initializing.")

    good_orbit(f, npts, Ttr = Ttr, max_attempts = max_attempts)
end

Main.ex-NormalisedPredictiveAsymmetryTest_logistic4.good_orbit

Generate time series

Okay, now we're ready to generate some time series.

# Example data
sys = logistic4
npts = 400

# Some initial conditons blow up, so iterate until we find a good one.
x, y, z, w = good_orbit(sys, npts, Ttr = 300)

p_ts = plot(xlabel = "Time step", ylabel = "Value", size = (382*2, 300),
    guidefont = font(13), tickfont = font(13), legendfont = font(13))
plot!(x, label = "x")
plot!(y, label = "y")
plot!(z, label = "z")
plot!(w, label = "w")

Test setup

Now, define a normalised predictive asymmetry test using a visitation frequency test to get our predictions for lags [-10:1; 1:10], symmetrically around the zero lag. To construct the partition of our state space reconstructions, we'll use the heuristic from ¹ that the number of intervals along each coordinate axis should not exceed $N^{\left(\frac{1}{D+1}\right)}$, where $N$ is the number of points in the time series and $D$ is the total dimension of the generalised reconstruction used for the transfer entropy analyses.

# Test setup
η_max = 10
k = 1; l = 1; m = 1 # embedding parameters, total dim = k + l + m
ηs = [-η_max:-1; 1:η_max]
bin = RectangularBinning(floor(Int, npts^(1/(k + l + m + 1))))
test = VisitationFrequencyTest(ηs = ηs, binning = bin,
    estimator = VisitationFrequency(b = 2)) # use base-2 logarithm
pa_test = NormalisedPredictiveAsymmetryTest(test, f = 1.0)

NormalisedPredictiveAsymmetryTest{VisitationFrequencyTest{20},10}(predictive_test = VisitationFrequencyTest{20}(k = 1, l = 1, m = 1, n = 1, τ = 1, estimator = VisitationFrequency(2), binning_summary_statistic = mean, binning = RectangularBinning(4), ηs = [-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]), f = 1.0)

That's all we need. Now, let's compute the predictive asymmetry statistic and see whether we find the expected causalities.

pas_xy = causality(x, y, pa_test)
pas_yx = causality(y, x, pa_test)
pas_yz = causality(y, z, pa_test)
pas_zy = causality(z, y, pa_test)
pas_zw = causality(z, w, pa_test)
pas_wz = causality(w, z, pa_test)
asymmetries = [pas_xy pas_yx pas_yz pas_zy pas_zw pas_wz]

10×6 Array{Float64,2}:
  0.982442   -1.5495   1.26956  -1.36158  0.402471   -0.7851 
  2.17342    -2.88352  2.46485  -2.18839  0.84751    -1.74353
  3.68144    -4.56551  3.55018  -3.13727  1.39781    -2.4449 
  5.04407    -6.00303  4.45423  -3.84393  1.91493    -3.84897
  6.12963    -7.22953  5.24629  -4.6367   2.54724    -4.9448 
  7.03143    -8.23725  5.75334  -4.88404  3.34744    -6.20234
  7.92002    -9.2731   6.56122  -5.4304   4.21761    -7.53951
  8.81343   -10.1881   7.35192  -5.87716  5.21972    -8.65381
  9.57043   -11.1987   8.14151  -6.07756  6.10343    -9.79629
 10.3148    -11.9987   8.80658  -6.4961   6.73666   -10.735  

Plot time series and results

Let's plot the results.

ymax = maximum(abs.(asymmetries)) * 1.05
p_pa = plot(xlabel = L"\eta \,\, [bits]", ylabel = L"\mathcal{A}(\eta) \,\, [bits]",
    ylims = (-ymax, ymax),  size = (382*2, 550), legend = :topleft,
    fg_legend = :transparent, bg_legend = :transparent,
    guidefont = font(13), tickfont = font(13), legendfont = font(12),
    xlims = (0.9, η_max+0.1))
hline!([1], ls = :dash, lc = :grey, label = "")
plot!(ηs[ηs .> 0], pas_xy, c = :black, label = L"x \to y")
plot!(ηs[ηs .> 0], pas_yz, c = :red, label = L"y \to z")
plot!(ηs[ηs .> 0], pas_zw, c = :blue, label = L"z \to w")
plot!(ηs[ηs .> 0], pas_yx, c = :black, ls = :dash, label = L"y \to x")
plot!(ηs[ηs .> 0], pas_zy, c = :red, ls = :dash, label = L"z \to y")
plot!(ηs[ηs .> 0], pas_wz, c = :blue, ls = :dash, label = L"w \to z")

Discussion

If the predictive asymmetry picks up on the causal interactions, it should be positive for the directions we know to be causal and negative in the directions that are known to be non-causal. Running this example for multiple initial conditions of the system, we find that this is very consistently the case!

For more information on about this method, wait patiently for our upcoming paper, which should be out as a preprint soon. Keep tuned for more!

---

1. Krakovská, A., Jakubík, J., Chvosteková, M., Coufal, D., Jajcay, N., & Paluš, M. (2018). Comparison of six methods for the detection of causality in a bivariate time series. Physical Review E, 97(4), 042207. https://journals.aps.org/pre/abstract/10.1103/PhysRevE.97.042207 ↩