Skip to content

Effect of dicretization scheme on transfer entropy estimates

Different ways of partitioning

The TransferOperatorGrid and VisitationFrequency transfer entropy estimators both operate on partitions on the delay reconstructions. Below, we demonstrate the four different ways of discretizing the state space.

First, let's load necessary packages and create some example time series.

using CausalityTools, Plots, Statistics

x = rand(300)
y = sin.(cumsum(rand(300)))*0.3 .+ 0.05*x*rand([-1.0, 1.0])

300-element Array{Float64,1}:
  0.13894516223376502 
  0.3218405110945099  
  0.3171438960238998  
  0.2984615391260007  
  0.1892228602740244  
  0.08071894552866779 
 -0.026999063158544098
 -0.155339709084643   
 -0.28162441015471495 
 -0.2691929013797143  
  ⋮                   
  0.12249576757296844 
  0.10542612932335023 
  0.07938771746540794 
  0.05510600690713506 
 -0.04168407832723132 
 -0.24313503812345189 
 -0.2893573283179723  
 -0.2199159013558739  
  0.006136882537549213

Now, create a generalised delay reconstruction. For computing the transfer entropy from x to y, we will create the generalised embedding vectors \{(y_{t+η}, y_{t}, y_{t-τ}, x_{t}) \}. In the opposite direction, we will have embedding vectors of the form \{(x_{t+η}, x_{t}, x_{t-τ}, y_{t})\}.

τ = 1 # embedding lag
η = 1 # forward prediction lag
E_xtoy = customembed(Dataset(x, y), Positions([2, 2, 2, 1]), Lags([η, 0, -τ, 0]))
E_ytox = customembed(Dataset(y, x), Positions([2, 2, 2, 1]), Lags([η, 0, -τ, 0]));

Now E_xtoy and E_ytox stores the lagged time series as a DynamicalSystems.Dataset, so that we can index the lagged variables by columns.

Okay, so we know which variables are in which columns and what lags they have. But the transfer entropy estimators don't have this information yet, so we need to provide it using a TEVars instance. This is necessary to organize the computation of marginal probabilities.

# Organize marginals
Tf = [1]     # future of target variable is in first column
Tpp = [2, 3] # present and past of target variable are in second and third columns
Spp = [4]    # the present of the source is in the fourth column.
v = TEVars(Tf, Tpp, Spp)

TEVars(Tf = [1], Tpp = [2, 3], Spp = [4], Cpp = Int64[])

Hyper-rectangles by subdivision of axes (ϵ::Int)

First, we use an integer number of subdivisions along each axis of the delay embedding when partitioning.

ϵs = 1:2:50 # integer number of subdivisions along each axis of the embedding
te_estimates_xtoy = zeros(length(ϵs))
te_estimates_ytox = zeros(length(ϵs))
vars = TEVars([1], [2, 3], [4])
estimator = VisitationFrequency(b = 2) # logarithm to base 2 gives units of bits

for (i, ϵ) in enumerate(ϵs)
    te_estimates_xtoy[i] = transferentropy(E_xtoy, vars, RectangularBinning(ϵ), estimator)
    te_estimates_ytox[i] = transferentropy(E_ytox, vars, RectangularBinning(ϵ), estimator)
end

p = plot()
plot!(p, ϵs, te_estimates_xtoy, label = "TE(x -> y)", lc = :black)
plot!(p, ϵs, te_estimates_ytox, label = "TE(y -> x)", lc = :red)
xlabel!(p, "# subdivisions along each axis")
ylabel!(p, "Transfer entropy (bits)")

Hyper-cubes of fixed size (ϵ::Float)

We do precisely the same, but use fixed-width hyper-cube bins. The values of the logistic map take values on [0, 1], so using bins width edge lengths 0.1 should give a covering corresponding to using 10 subdivisions along each axis of the delay embedding. We let ϵ take values on [0.05, 0.5].

ϵs = 0.02:0.02:0.5
te_estimates_xtoy = zeros(length(ϵs))
te_estimates_ytox = zeros(length(ϵs))
vars = TEVars([1], [2, 3], [4])
estimator = VisitationFrequency(b = 2)

for (i, ϵ) in enumerate(ϵs)
    te_estimates_xtoy[i] = transferentropy(E_xtoy, vars, RectangularBinning(ϵ), estimator)
    te_estimates_ytox[i] = transferentropy(E_ytox, vars, RectangularBinning(ϵ), estimator)
end

plot(ϵs, te_estimates_xtoy, label = "TE(x -> y)", lc = :black)
plot!(ϵs, te_estimates_ytox, label = "TE(y -> x)", lc = :red)
xlabel!("Hypercube edge length")
ylabel!("Transfer entropy (bits)")
xflip!()

Hyper-rectangles of fixed size (ϵ::Vector{Float})

It is also possible to use hyper-rectangles, by specifying the edge lengths along each coordinate axis of the delay embedding. In our case, we use a four-dimensional, embedding, so we must provide a 4-element vector of edge lengths

# Define slightly different edge lengths along each axis
ϵs_x1 = LinRange(0.05, 0.5, 10)
ϵs_x2 = LinRange(0.02, 0.4, 10)
ϵs_x3 = LinRange(0.08, 0.6, 10)
ϵs_x4 = LinRange(0.10, 0.3, 10)

te_estimates_xtoy = zeros(length(ϵs_x1))
te_estimates_ytox = zeros(length(ϵs_x1))
vars = TEVars([1], [2, 3], [4])
estimator = VisitationFrequency(b = 2)

mean_ϵs = zeros(10)

for i ∈ 1:10
    ϵ = [ϵs_x1[i], ϵs_x2[i], ϵs_x3[i], ϵs_x4[i]]
    te_estimates_xtoy[i] = transferentropy(E_xtoy, vars, RectangularBinning(ϵ), estimator)
    te_estimates_ytox[i] = transferentropy(E_ytox, vars, RectangularBinning(ϵ), estimator)

    # Store average edge length (for plotting)
    mean_ϵs[i] = mean(ϵ)
end

plot(mean_ϵs, te_estimates_xtoy, label = "TE(x -> y)", lc = :black)
plot!(mean_ϵs, te_estimates_ytox, label = "TE(y -> x)", lc = :red)
xlabel!("Average hypercube edge length")
ylabel!("Transfer entropy (bits)")
xflip!()

Hyper-rectangles by variable-width subdivision of axes (ϵ::Vector{Int})

Another way to construct hyper-rectangles is to subdivide each coordinate axis into segments of equal length. In our case, we use a four-dimensional, embedding, so we must provide a 4-element vector providing the number of subdivisions we want along each axis.

# Define different number of subdivisions along each axis.
ϵs = 3:50
mean_ϵs = zeros(length(ϵs))
estimator = VisitationFrequency(b = 2)

te_estimates_xtoy = zeros(length(ϵs))
te_estimates_ytox = zeros(length(ϵs))
vars = TEVars([1], [2, 3], [4])

for (i, ϵᵢ) ∈ enumerate(ϵs)
    ϵ = [ϵᵢ - 1, ϵᵢ, ϵᵢ, ϵᵢ + 1]
    te_estimates_xtoy[i] = transferentropy(E_xtoy, vars, RectangularBinning(ϵ), estimator)
    te_estimates_ytox[i] = transferentropy(E_ytox, vars, RectangularBinning(ϵ), estimator)

    # Store average number of subdivisions for plotting
    mean_ϵs[i] = mean(ϵ)
end

plot(mean_ϵs, te_estimates_xtoy, label = "TE(x -> y)", lc = :black)
plot!(mean_ϵs, te_estimates_ytox, label = "TE(y -> x)", lc = :red)
xlabel!("Average number of subdivisions along the embedding axes")
ylabel!("Transfer entropy (bits)")