Causality in dynamical systems
General syntax
The syntax for quantifying directional causality between variables of dynamical systems is
causality(system::DynamicalSystem, setup::DynamicalSystemSetup, test::CausalityTest)
Here,
systemis either aDiscreteDynamicalSystemorContinuousDynamicalSysteminstance. SeeDynamicalSystems.jlfor more info on how to construct and initialise dynamical systems,setupis an instance of eitherContinuousSystemSetuporDiscreteSystemSetup. This gives the information on trajectory lengths, sampling steps and solver options. It also instructs the causality estimators about which variables of the dynamical system from which the time series forsourceandtargetare to be sampled,testis an instance of a causality test.
What happens under the hood is that the dynamical system you provided is iterated using the provided setup parameters. A trajectory for the system is then recorded and subsampled according to your specifications in setup. Then, the causality test is applied to the time series for the variables you indicated in the setup.
Continuous systems
#
CausalityToolsBase.causality — Method.
causality(system::ContinuousDynamicalSystem, setup::ContinuousSystemSetup,
test::CausalityTest)
Apply the causality test to the given continuous system using the provided setup parameters.
Example
Let's re-create the lorenz_lorenz system that ships with CausalityTools and analyse the dynamical influence between some of its components.
function eom_lorenz_lorenz(u, p, t)
c_xy, c_yx, a₁, a₂, a₃, b₁, b₂, b₃ = (p...,)
x1, x2, x3, y1, y2, y3 = (u...,)
dx1 = -a₁*(x1 - x2) + c_yx*(y1 - x1)
dx2 = -x1*x3 + a₂*x1 - x2
dx3 = x1*x2 - a₃*x3
dy1 = -b₁*(y1 - y2) + c_xy*(x1 - y1)
dy2 = -y1*y3 + b₂*y1 - y2
dy3 = y1*y2 - b₃*y3
return SVector{6}(dx1, dx2, dx3, dy1, dy2, dy3)
end
function lorenz_lorenz(; u0 = rand(6),
c_xy = 0.2, c_yx = 0.0,
a₁ = 10, a₂ = 28, a₃ = 8/3,
b₁ = 10, b₂ = 28, b₃ = 9/3)
ContinuousDynamicalSystem(eom_lorenz_lorenz, u0, [c_xy, c_yx, a₁, a₂, a₃, b₁, b₂, b₃])
end
# Create an instance of the system with default parameters
sys = lorenz_lorenz()
# Analysis setup. We'll check for an influence from variable x1 to y1, which are the
# 1st and 4st components of the orbit. We'll generate orbits consisting of 500 points.
setup = ContinuousSystemSetup(source = 1, target = 4, n_pts = 500)
# Predictive asymmetry causality test
predtest = TransferOperatorGridTest(binning = RectangularBinning(5), ηs = -5:5)
test_pa = PredictiveAsymmetryTest(predictive_test = predtest)
causality(sys, setup, test_pa)
Discrete systems
#
CausalityToolsBase.causality — Method.
causality(system::DiscreteDynamicalSystem, setup::DiscreteSystemSetup,
test::CausalityTest)
Apply the causality test to the given discrete system using the provided setup parameters.
Example
# Use one of the built-in dynamical systems that ship with CausalityTools,
# `logistic2_unidir`, which is a system of two chaotic logistic maps `x` and `y` with
# coupling `x` to `y`, with `c_xy` controlling the coupling strength.
sys = logistic2_unidir(c_xy = 1.0)
# We treat the first variable as the source and the second variable as the target.
# Iterate the system until we have time series with 200 observations.
setup = DiscreteSystemSetup(source = 1, target = 2, n_pts = 200)
# Define some causality tests (with default values, you'd want to tweak the
# parameters to your needs)
test_ccm = ConvergentCrossMappingTest(timeseries_lengths = [45, 50], n_reps = 20)
test_cm = CrossMappingTest(n_reps = 10)
test_vf = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = 1:5)
test_tog = TransferOperatorGridTest(binning = RectangularBinning(5), ηs = 1:5)
test_jdd = JointDistanceDistributionTest()
test_jddt = JointDistanceDistributionTTest()
predtest = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = -5:5)
test_pa = PredictiveAsymmetryTest(predictive_test = predtest)
# Compute the causality statistics using the different tests
causality(sys, setup, test_ccm)
causality(sys, setup, test_cm)
causality(sys, setup, test_vf)
causality(sys, setup, test_tog)
causality(sys, setup, test_jdd)
causality(sys, setup, test_jddt)
Setting up the analysis
The setup for continuous and discrete systems is slightly different, so there are two separate types where you can give instructions about step length, time series length, etc, to the solvers.
Continuous system setup
#
CausalityTools.IntegrationDynamicalSystems.ContinuousSystemSetup — Type.
ContinuousSystemSetup(resampling = NoResampling, dt::Int = 0.01, Ttr::Int = 0,
sample_step::Int = 1, n_pts::Int = 100,
diffeq = (alg = SimpleDiffEq.SimpleATsit5(), abstol = 1.0e-6, reltol = 1.0e-6),
source, target)
Setup for causality analysis of continuous dynamical systems.
Mandatory keyword arguments
source: The variable(s) to use assourcewhen callingcausality. Usually integer(s).target: The variable(s) to use astargetwhen callingcausality. Usually integer(s).
Optional keywords
resampling: An instance of a resampling scheme. Defaults toNoResampling().dt::Int: The time step when iterating the system. Defaults to 0.01.Ttr::Int: The number of transient iterations before starting sampling. Defaults to 0.sample_step::Int: The sampling step in the final orbit. Ifsample_step> 1, then the orbit is generated until timeT*sample_step, after which eachsample_stepth point is drawn to generate the final orbit.diffeq: Arguments propagated to init DifferentialEquations.jl. Defaults todiffeq = (alg = SimpleDiffEq.SimpleATsit5(), abstol = 1.0e-6, reltol = 1.0e-6).n_pts::Int: The number of points in the orbit. Defaults to 100.
Discrete system setup
#
CausalityTools.IntegrationDynamicalSystems.DiscreteSystemSetup — Type.
DiscreteSystemSetup(resampling = NoResampling, dt::Int = 1, Ttr::Int = 0,
diffeq = (alg = SimpleDiffEq.SimpleATsit5(), abstol = 1.0e-6, reltol = 1.0e-6),
n_pts = 100, source, target)
Setup for causality analysis of discrete dynamical systems.
Mandatory keyword arguments
source: The variable(s) to use assourcewhen callingcausality. Usually integer(s).target: The variable(s) to use astargetwhen callingcausality. Usually integer(s).
Optional keywords
resampling: An instance of a resampling scheme. Defaults toNoResampling().dt::Int: The time step when iterating the system.Ttr::Int: The number of transient iterations before starting sampling.diffeq: Arguments propagated to init DifferentialEquations.jl. Defaults todiffeq = (alg = SimpleDiffEq.SimpleATsit5(), abstol = 1.0e-6, reltol = 1.0e-6).n_pts::Int: The number of points in the orbit.