Predictive asymmetry

CausalityTools.PredictiveAsymmetry.predictive_asymmetryFunction

General interface

predictive_asymmetry(s, t, [c], 
    estimator::TransferEntropyEstimator, ηs; 
    d𝒯 = 1, dT = 1, dS = 1, τT = -1, τS = -1, 
    [dC = 1, τC = -1,],
    normalize::Bool = false, f::Real = 1.0) → Vector{Float64}

Compute the predictive asymmetry[Haaga2020] 𝔸(st) for source time series s and target time series t over prediction lags ηs, using the given estimator and embedding parameters d𝒯, dT, dS, τT, τS.

If a conditional time series c is provided, compute 𝔸(stc). Then, dC and τC controls the embedding dimension and embedding lag for the conditional variable.

Returns

Returns a vector containing the predictive asymmetry for each value of ηs.

Normalization (hypothesis test)

If normalize == true (the default), then compute the normalized predictive asymmetry 𝒜.

In this case, for each $\eta$ in ηs, compute 𝒜(η) by normalizing 𝔸(η) to some fraction f of the mean transfer entropy over prediction lags $-\eta, ..., \eta$ (exluding lag 0).

Haaga et al. (2020)[Haaga2020] uses a normalization with f=1.0 as a built-in hypothesis test, avoiding more computationally costly surrogate testing.

Estimators

Any estimator that works for transferentropy will also work with predictive_asymmetry. It is recommended to use either the rectangular binning-based methods or the symbolic estimators for the fastest computations.

Examples

using CausalityTools 

# Some example time series
x, y = rand(100), rand(100)

# 𝔸(x → y) over prediction lags 1:5
𝔸reg  = predictive_asymmetry(x, y, VisitationFrequency(RectangularBinning(5)), 1:5) 
source

Example: no coupling

Here, we'll compute the predictive asymmetry on 100 different sets of random time series. Because there is no dynamical coupling between the time series, we expect the predictive asymmetry to be zero.

We'll use a visitation frequency estimator.

using CausalityTools, Plots, StatsBase

# Define prediction lags and estimator
ηs = 1:10
est = VisitationFrequency(RectangularBinning(5)) # guess that 5 bins along each coordinate axis is sufficient

nreps = 100
𝔸xy = zeros(nreps, length(ηs))
𝔸yx = zeros(nreps, length(ηs))

for i = 1:nreps
    # Some example time series
    x, y = rand(1000), rand(1000)
    𝔸xy[i, :] = predictive_asymmetry(x, y, est, ηs)
    𝔸yx[i, :] = predictive_asymmetry(y, x, est, ηs)
end

plot()
plot!(ηs, dropdims(mean(𝔸xy, dims = 1), dims = 1), label = "X to Y", c = :black)
plot!(ηs, dropdims(mean(𝔸yx, dims = 1), dims = 1), label = "Y to X", c = :red)
ylims!((-0.1, 0.1))
xlabel!("Prediction lag (η)")
ylabel!("Predictive asymmetry (bits)")

As expected, because there is no dynamical coupling between the variables, the predictive asymmetry is around zero for all prediction lags.

Example: coupled Ulam maps

Here, we'll use the same example as in Schreiber's transfer entropy paper and see if the predictive asymmetry is able to decipher the directional coupling in the Ulam system.

Let's start by defining the Ulam system. Then, we generate time series for different values of the parameter ε and compute the predictive asymmetry for prediction lags ηs = 1:10.

We expect that the predictive asymmetry is positive if there is a directional coupling, and zero or negative if there is no directional coupling.

using DynamicalSystems, CausalityTools, Plots, Random, StatsBase, TimeseriesSurrogates

Random.seed!(12234)

function ulam_system(dx, x, p, t)
    f(x) = 2 - x^2
    ε = p[1]
    dx[1] = f(ε*x[length(dx)] + (1-ε)*x[1])
    for i in 2:length(dx)
        dx[i] = f(ε*x[i-1] + (1-ε)*x[i])
    end
end

ds = DiscreteDynamicalSystem(ulam_system, rand(100) .- 0.5, [0.04])
trajectory(ds, 1000; Ttr = 1000)

εs = 0.02:0.02:1.0
ηs = 1:10
pas_x1x2 = zeros(length(εs), length(ηs));
pas_x2x1 = zeros(length(εs), length(ηs))

for (i, ε) in enumerate(εs)
    set_parameter!(ds, 1, ε)
    # Use time series consisting of 1000 points
    tr = trajectory(ds, 1000; Ttr = 5000)
    X1 = tr[:, 1]; X2 = tr[:, 2]
    @assert !any(isnan, X1)
    @assert !any(isnan, X2)
    binning = RectangularBinning(0.2) # guess an appropriate bin width of 0.2
    pas_x1x2[i, :] = predictive_asymmetry(X1, X2, VisitationFrequency(binning), ηs)
    pas_x2x1[i, :] = predictive_asymmetry(X2, X1, VisitationFrequency(binning), ηs)
end

Let's pick the predictive asymmetry at prediction lag η = 10.

# Pick predictive asymmetry at η = 10
pa_x1x2 = [pas_x1x2[i, end] for i = 1:length(εs)]
pa_x2x1 = [pas_x2x1[i, end] for i = 1:length(εs)]

plot()
plot(εs, pa_x1x2, label = "X1 to X2", c = :black, lw = 1.5)
plot!(εs, pa_x2x1, label = "X2 to X1", c = :red)
xlabel!("epsilon")
ylabel!("Predictive asymmetry (bits)")
hline!([0], ls = :dash, label = "", c = :grey)

As expected, predictive asymmetry is positive in the direction $X_1 \\to X_2$, and negative in the direction $X_2 \\to X_1$.

Are the results significant? We can test that using a surrogate test. We'll use random shuffle surrogates from the TimeseriesSurrogates.jl package.

Using the same parameters, as before, we generate time series for the same values of εs, but for each ε, we generate 100 different shuffled versions of the source time series. For each of those surrogate time series, we compute the predictive asymmetry and pick the value at η = 10

n_surr = 100
pas_x1x2_surr = zeros(length(εs), n_surr);
pas_x2x1_surr = zeros(length(εs), n_surr)

for (i, ε) in enumerate(εs)
    set_parameter!(ds, 1, ε)
    tr = trajectory(ds, 1000; Ttr = 5000)
    X1 = tr[:, 1]; @assert !any(isnan, X1)
    X2 = tr[:, 2]; @assert !any(isnan, X2)
    S1 = surrogenerator(X1, RandomShuffle())
    S2 = surrogenerator(X2, RandomShuffle())

    binning = RectangularBinning(0.2) # guess an appropriate bin width of 0.2

    pas_x1x2_surr[i, :] = [predictive_asymmetry(S1(), X2, VisitationFrequency(binning), ηs)[end] for k = 1:n_surr]
    pas_x2x1_surr[i, :] = [predictive_asymmetry(S2(), X1, VisitationFrequency(binning), ηs)[end] for k = 1:n_surr]
end

We'll plot the 95th percentile of the surrogate ensemble predictive asymmetries (dashed lines) along with the original predictive asymmetries (solid lines). If the predictive asymmetry computed on the original time series exceed that of the 95th percentile of the surrogate ensemble, then the predictive asymmetry is significant.

uq_x1x2 = [quantile(pas_x1x2_surr[i, :], 0.95) for i = 1:length(εs)]
uq_x2x1 = [quantile(pas_x2x1_surr[i, :], 0.95) for i = 1:length(εs)]
ymax = maximum(abs.([pa_x1x2; pa_x2x1; uq_x1x2; uq_x2x1]))*1.1
plot(ylims = (-ymax, ymax))
plot!(εs, pa_x1x2, label = "X1 to X2", c = :black, lw = 1.5)
plot!(εs, uq_x1x2, label = "X1 to X2", c = :black, ls = :dot)
plot!(εs, pa_x2x1, label = "X2 to X1", c = :red)
plot!(εs, uq_x2x1, label = "X2 to X1", c = :red, ls = :dot)

xlabel!("epsilon")
ylabel!("Predictive asymmetry (bits)")
hline!([0], ls = :dash, label = "", c = :grey)

The predictive asymmetry in the direction $X_1 \\to X_2$ is above the 95th percentile of the surrogates for most parameters. Thus, we reject the null hypothesis that $X_1$ does not influence $X_2$ and accept the alternative hypothesis that $X_1$ does in fact influence $X_2$. In the context of causal inference, this means that there is evidence of directional dynamical coupling from $X_1$ to $X_2$ and that $X_1$ is causally coupled to $X_2$.

In the oppposite direction, $X_2 \\to X_1$, where there should not be any directional coupling (by definition), the predictive asymmetries are below the surrogates, so we cannot reject the null hypothesis that $X_2$ does not influence $X_1$. Hence, there is no evidence of directional dynamical coupling from $X_2$ to $X_1$.

  • Haaga2020Haaga, Kristian Agasøster, David Diego, Jo Brendryen, and Bjarte Hannisdal. "A simple test for causality in complex systems." arXiv preprint arXiv:2005.01860 (2020).