

SMeasureTest

# CausalityTools.CausalityTests.SMeasureTest — Type.

SMeasureTest(m, τ, K)
SMeasureTest(m::Int, τ, K; metric::Metric = SqEuclidean(), 
    tree_metric::Metric = Euclidean())

S-measure test [1] for the directional dependence between data series.

Algorithm

1. Create an m-dimensional embeddings of both x and y, resulting in N different m-dimensional embedding points $X = \{x_1, x_2, \ldots, x_N \}$ and $X = \{y_1, y_2, \ldots, y_N \}$. τ controls the embedding lag.
2. Let $r_{i,j}$ and $s_{i,j}$ be the indices of the k-th nearest neighbors of $x_i$ and $y_i$, respectively.
3. Compute the the mean squared Euclidean distance to the $k$ nearest neighbors for each $x_i$, using the indices $r_{i, j}$.

$$R_i^{(k)}(x) = \dfrac{1}{k} \sum_{i=1}^{k}(x_i, x_{r_{i,j}})^2$$

• Compute the y-conditioned mean squared Euclidean distance to the $k$ nearest neighbors for each $x_i$, now using the indices $s_{i,j}$.

$$R_i^{(k)}(x|y) = \dfrac{1}{k} \sum_{i=1}^{k}(x_i, x_{s_{i,j}})^2$$

• Define the following measure of independence, where $0 \leq S \leq 1$, and low values indicate independence and values close to one occur for synchronized signals.

$$S^{(k)}(x|y) = \dfrac{1}{N} \sum_{i=1}^{N} \dfrac{R_i^{(k)}(x)}{R_i^{(k)}(x|y)}$$

Examples

# Initialise test, specifying embedding dimension, embedding lag 
# and number of nearest neighbors
test = SMeasureTest(m = 4, τ = 3, K = 2:10)

References

1. Quian Quiroga, R., Arnhold, J. & Grassberger, P. [2000] “Learning driver-response relationships from synchronization patterns,” Phys. Rev. E61(5), 5142–5148.

source

Example

Time series

First, create an orbit of the built-in unidirectionally coupled henon2 map system, and a set of random time series for comparison.

using CausalityTools, DynamicalSystems, Plots, Distributions

npts, Ttr = 5000, 500
x, y = columns(trajectory(henon2(c_xy = 1.0), npts - 1, Ttr = Ttr))
xr, yr = rand(Uniform(-1, 1), npts), rand(Uniform(-1, 1), npts)
[x y xr yr]

5000×4 Array{Float64,2}:
  0.394001  -1.21518    0.305167      0.692067  
  0.880209   0.394001   0.762157     -0.442792  
  0.743432   0.880209  -0.242306      0.477283  
  1.11137    0.743432  -0.591254      0.766874  
  0.387882   1.11137    0.0867659    -0.0211271 
  1.58296    0.387882  -0.000578745   0.573234  
 -0.989394   1.58296    0.393555     -0.0184107 
  0.895988  -0.989394   0.644794     -0.837293  
  0.300387   0.895988   0.54911       0.203863  
  1.57856    0.300387   0.12882       0.801228  
  ⋮                                             
  1.5317     0.39208    0.249539      0.705748  
 -0.828472   1.5317     0.261045     -0.196222  
  1.17314   -0.828472   0.457707      0.0560364 
 -0.224807   1.17314   -0.412061     -0.485966  
  1.7014    -0.224807  -0.583629     -0.517048  
 -1.56222    1.7014    -0.0360618     0.00988211
 -0.530112  -1.56222    0.894293     -0.900317  
  0.650315  -0.530112  -0.46669      -0.645396  
  0.818057   0.650315  -0.0960493    -0.750436  

Let's plot the first 100 points of each time series.

p_det = plot(xlabel = "", ylabel = "Value", title = "Coupled Henon maps")
plot!(x[1:100], label = "x", marker = stroke(:black), c = :black)
plot!(y[1:100], label = "y", marker = stroke(:red), c = :red)
p_rand = plot(xlabel = "Time step", ylabel = "Value", title = "Random time series")
plot!(xr[1:100], label = "xr", c = :blue)
plot!(yr[1:100], label = "yr", c = :purple, ls = :dash)

plot(p_det, p_rand, layout = grid(2, 1), size = (382*2, 400), legend = :bottomright,
    tickfont = font(13), guidefont = font(13), legendfont = font(13))

Test setup

Create a SMeasureTest, which contains the test parameters. These are the embedding dimension m, the embedding lag τ, and the number of nearest neighbors K. We'll compute the test with fixed embedding, but a varying number of nearest neighbors ks = 2:10.

ks = 2:10
test = SMeasureTest(m = 4, τ = 1, K = ks)

SMeasureTest{9,Distances.SqEuclidean,Euclidean}(m = 4, τ = 1, K = 2:10, metric = Distances.SqEuclidean(0.0), tree_metric = Euclidean(0.0))

Analysis

Compute S-measure statistic separately in both directions, both for the random time series, and for the Henon maps. The test will return a vector of length ks.

Ss_r_xy = causality(xr, yr, test)
Ss_r_yx = causality(yr, xr, test)
Ss_henon_xy = causality(x, y, test)
Ss_henon_yx = causality(y, x, test);
[Ss_r_xy Ss_r_yx Ss_henon_xy Ss_henon_yx]

9×4 Array{Float64,2}:
 0.0158544  0.0160419  0.817326  0.79036 
 0.0144574  0.0145955  0.728242  0.667714
 0.0158563  0.0158216  0.689484  0.614211
 0.017319   0.017317   0.671805  0.58075 
 0.0188464  0.0187807  0.655392  0.557013
 0.0202119  0.0201693  0.646425  0.548321
 0.021545   0.0215229  0.641744  0.541835
 0.0228955  0.0227664  0.638436  0.539838
 0.0240551  0.0239909  0.636126  0.533974

The first two columns contain the s-measures computed for our deterministic Henon map time series, and the two last columns contain the s-measures computed for the stochastic time series.

Results

Let's plot the results.

plot(xlabel = "# nearest neighbors (k)", ylabel = "S", ylims = (-0.05, 1.05))
plot!(ks, Ss_r_xy,  label = "random uncoupled system (x -> y)", marker = stroke(2), c = :black)
plot!(ks, Ss_r_yx,  label = "random uncoupled system (y -> x)", marker = stroke(2), c = :red)
plot!(ks, Ss_henon_xy, marker = stroke(2), label = "henon unidir (x -> y)")
plot!(ks, Ss_henon_yx, marker = stroke(2), label = "henon unidir (y -> x)")

Discussion

For uncoupled time series, we expect the value of $S$ to be close to zero. For strongly coupled time series, the value of $S$ should be nonzero and approaching 1. This is exactly what we get: for time random time series, the value of $S$ is close to zero and for the Henon maps, it's clearly non-zero.