Tutorial: SMeasureTest
This tutorials gives an example of how to use the SMeasureTest.
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}:
1.70356 0.0716883 -0.7584 0.247085
-1.4806 1.70356 -0.708043 0.948896
-0.281112 -1.4806 -0.0256723 -0.372124
0.876796 -0.281112 -0.688912 -0.922445
0.546895 0.876796 -0.799998 0.0392815
1.36394 0.546895 -0.551733 -0.088644
-0.296276 1.36394 0.3714 -0.54769
1.7214 -0.296276 0.704878 0.500032
-1.65211 1.7214 0.924478 0.31865
-0.813061 -1.65211 0.512393 0.586339
⋮
0.474355 -1.19028 0.679082 0.924796
0.817905 0.474355 0.951915 0.725182
0.873339 0.817905 -0.0188547 0.589741
0.882651 0.873339 -0.70909 -0.841099
0.882929 0.882651 0.868487 0.316532
0.885233 0.882929 -0.90695 0.900742
0.881242 0.885233 0.676192 0.282829
0.888982 0.881242 0.143864 -0.706657
0.874083 0.888982 -0.460734 -0.726076
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",a ylabel = "Value", title = "Random time series")
plot!(xr[1:100], label = "xr", c = :blue)
plot!(yr[1:100], label = "yr", c = :purple)
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.0165754 0.0161279 0.794415 0.770415
0.0147225 0.0145271 0.710664 0.656022
0.0158688 0.0156643 0.675702 0.587412
0.0172917 0.0171041 0.650984 0.561205
0.0187841 0.0186656 0.636835 0.540516
0.020161 0.020076 0.627252 0.538196
0.0214621 0.0214282 0.624118 0.532514
0.0227114 0.0227224 0.62149 0.528903
0.0238806 0.0239897 0.620613 0.5255
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.