21  Counterfactual Dynamics: Alternative Trajectories

Status: Draft

v0.4

21.1 Learning Objectives

After reading this chapter, you will be able to:

  • Define counterfactual trajectories in dynamical systems
  • Compute unit-level counterfactual dynamics using shared exogenous noise
  • Understand bounds and partial identification for dynamic counterfactuals
  • Apply counterfactual reasoning to dynamical systems

21.2 Introduction

This chapter explores counterfactual reasoning in dynamical systems—what alternative trajectories would have been possible for specific units under different interventions, given the same exogenous noise realisation. This is part of Imagining in the Dynamical world—what alternative dynamic trajectories are possible?

21.3 Counterfactual Trajectories in Dynamical Systems

21.3.1 Definition

For a fixed exogenous realisation \(\mathbf{u}\), the counterfactual trajectory under intervention \(\iota\) in a dynamical system is:

\[ Y^{\iota}_{1:T}(\mathbf{u}) \]

Interpretation: “Same unit (same \(\mathbf{u}\)), different world (different intervention), same dynamic system.”

21.3.2 Key Requirement

Counterfactuals in dynamical systems require shared exogenous noise: - Same unit \(\Leftrightarrow\) same \(\mathbf{u}\) (exogenous noise realisation) - Different worlds \(\Leftrightarrow\) different interventions - Same dynamics \(\Leftrightarrow\) same system structure

21.4 Computing Counterfactual Dynamics

21.4.1 Unit-Level Approach

Steps: 1. Infer \(\mathbf{u}\): From observed trajectory, infer the exogenous noise realisation for this unit 2. Simulate counterfactual: With same \(\mathbf{u}\), simulate under alternative intervention 3. Compare trajectories: Compare counterfactual trajectory to observed trajectory

Challenge: Inferring \(\mathbf{u}\) from observations requires strong assumptions (full structural model with known noise structure).

21.4.2 Implementation: Computing Counterfactual Trajectories

Here’s an example showing how to compute counterfactual trajectories for a specific unit:

# Find project root and include ensure_packages.jl
project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(project_root, "scripts", "ensure_packages.jl"))

@auto_using DifferentialEquations Random CairoMakie

Random.seed!(123)

# Example: Treatment timing counterfactual
# Observed: Treatment started at t=10
# Counterfactual: What if treatment started at t=5?

β = 0.2  # Disease progression
α = 0.3  # Treatment effectiveness

# Step 1: Simulate observed trajectory (treatment at t=10)
function disease_observed!(du, u, p, t)
    """Disease model with observed treatment trajectory (treatment at t=10)."""
    X = u[1]
    # Treatment starts at t=10
    A = t >= 10.0 ? 1.0 : 0.0
    du[1] = β * X - α * A * X
end

u0_observed = [0.1]  # Initial severity
tspan = (0.0, 20.0)
prob_observed = ODEProblem(disease_observed!, u0_observed, tspan)
sol_observed = solve(prob_observed, Tsit5())

# Step 2: Infer the "unit" (in this deterministic case, unit = initial condition)
# For deterministic systems, the unit is fully determined by initial conditions
# For stochastic systems, we'd need to infer the noise realisation

# Step 3: Simulate counterfactual (same unit, different intervention: treatment at t=5)
function disease_counterfactual!(du, u, p, t)
    """Disease model with counterfactual treatment trajectory (treatment at t=5)."""
    X = u[1]
    # Counterfactual: Treatment starts at t=5 (earlier)
    A = t >= 5.0 ? 1.0 : 0.0
    du[1] = β * X - α * A * X
end

# Same initial condition (same unit)
u0_counterfactual = u0_observed
prob_counterfactual = ODEProblem(disease_counterfactual!, u0_counterfactual, tspan)
sol_counterfactual = solve(prob_counterfactual, Tsit5())

# Visualise comparison
let
fig = Figure(size = (800, 400))
ax = Axis(fig[1, 1], title = "Counterfactual Dynamics: Treatment Timing", 
          xlabel = "Time", ylabel = "Disease Severity")

lines!(ax, sol_observed.t, [u[1] for u in sol_observed.u], 
       label = "Observed (treatment at t=10)", linewidth = 2, color = :blue)
lines!(ax, sol_counterfactual.t, [u[1] for u in sol_counterfactual.u], 
       label = "Counterfactual (treatment at t=5)", linewidth = 2, color = :red, linestyle = :dash)

# Mark intervention times
vlines!(ax, [5.0, 10.0], color = :gray, linestyle = :dot, linewidth = 1)
axislegend(ax)

    fig  # Only this gets displayed
end

# Compare outcomes
final_observed = sol_observed.u[end][1]
final_counterfactual = sol_counterfactual.u[end][1]
println("Final severity:")
println("  Observed (treatment at t=10): ", round(final_observed, digits=3))
println("  Counterfactual (treatment at t=5): ", round(final_counterfactual, digits=3))
println("  Improvement from earlier treatment: ", round((final_observed - final_counterfactual) / final_observed * 100, digits=1), "%")
Final severity:
  Observed (treatment at t=10): 0.264
  Counterfactual (treatment at t=5): 0.064
  Improvement from earlier treatment: 75.9%

21.4.3 Example: “What Would Have Happened If We Had Intervened Earlier?”

Question: For a specific patient, what would have happened if we had started treatment at time \(t_0\) instead of \(t_1\)?

Process: 1. Infer \(\mathbf{u}\) from observed trajectory (treatment started at \(t_1\)) 2. Simulate counterfactual: same \(\mathbf{u}\), but treatment starts at \(t_0\) 3. Compare: counterfactual trajectory vs observed trajectory

21.5 Bounds and Partial Identification

When full counterfactual dynamics are not identified, we can still obtain bounds:

  • Non-parametric bounds: Range of possible counterfactual trajectories
  • Sensitivity parameters: How results change with assumptions about unobserved confounders
  • Uncertainty propagation: How uncertainty in \(\mathbf{u}\) affects counterfactual trajectories

21.5.1 Implementation: Bounds for Counterfactual Dynamics

When we cannot fully identify counterfactual trajectories, we can compute bounds:

# Find project root and include ensure_packages.jl
project_root = let
    current = pwd()
    while !isfile(joinpath(current, "Project.toml")) && !isfile(joinpath(current, "_quarto.yml"))
        parent = dirname(current)
        parent == current && break
        current = parent
    end
    current
end
include(joinpath(project_root, "scripts", "ensure_packages.jl"))

@auto_using DifferentialEquations CairoMakie

# Example: Counterfactual with uncertainty in treatment effect
# We know treatment helps, but exact effect α is uncertain: α ∈ [0.2, 0.4]

β = 0.2  # Disease progression (known)

# Lower bound: α = 0.2 (weakest treatment effect)
function disease_lower!(du, u, p, t)
    """Disease model: lower bound (α = 0.2) for partial identification."""
    X = u[1]
    A = t >= 5.0 ? 1.0 : 0.0
    du[1] = β * X - 0.2 * A * X  # Lower bound on treatment effect
end

# Upper bound: α = 0.4 (strongest treatment effect)
function disease_upper!(du, u, p, t)
    """Disease model: upper bound (α = 0.4) for partial identification."""
    X = u[1]
    A = t >= 5.0 ? 1.0 : 0.0
    du[1] = β * X - 0.4 * A * X  # Upper bound on treatment effect
end

u0 = [0.1]
tspan = (0.0, 20.0)

prob_lower = ODEProblem(disease_lower!, u0, tspan)
sol_lower = solve(prob_lower, Tsit5())

prob_upper = ODEProblem(disease_upper!, u0, tspan)
sol_upper = solve(prob_upper, Tsit5())

# Visualise bounds
let
fig = Figure(size = (800, 400))
ax = Axis(fig[1, 1], title = "Counterfactual Bounds: Treatment Effect Uncertainty", 
          xlabel = "Time", ylabel = "Disease Severity")

# Shade the region between bounds
# Solutions may have different numbers of time points due to adaptive stepping.
# Use the minimum length to ensure arrays align for plotting.
n_points = min(length(sol_lower.t), length(sol_upper.t))
t_common = sol_lower.t[1:n_points]
y_lower = [sol_lower.u[i][1] for i in 1:n_points]
y_upper = [sol_upper.u[i][1] for i in 1:n_points]
band!(ax, t_common, y_lower, y_upper, color = (:red, 0.3))

lines!(ax, t_common, y_lower, 
       label = "Lower bound (α=0.2)", linewidth = 2, color = :red, linestyle = :dash)
lines!(ax, t_common, y_upper, 
       label = "Upper bound (α=0.4)", linewidth = 2, color = :red, linestyle = :dash)

axislegend(ax)

    fig  # Only this gets displayed
end

println("Counterfactual bounds (final severity):")
println("  Lower bound (weakest treatment): ", round(sol_lower.u[end][1], digits=3))
println("  Upper bound (strongest treatment): ", round(sol_upper.u[end][1], digits=3))
println("  Range: ", round(sol_upper.u[end][1] - sol_lower.u[end][1], digits=3))
Counterfactual bounds (final severity):
  Lower bound (weakest treatment): 0.26
  Upper bound (strongest treatment): 0.013
  Range: -0.247

21.6 World Context

This chapter addresses Imagining in the Dynamical world—what alternative dynamic trajectories are possible? Counterfactual dynamics explore alternative trajectories in dynamical systems—what would have happened for specific units under different interventions, given the same exogenous noise realisation.

21.7 Key Takeaways

  1. Counterfactual trajectories: Alternative dynamic trajectories for specific units
  2. Shared exogenous noise: Maintains unit identity across counterfactual worlds
  3. Unit-level reasoning: Counterfactuals are for specific units, not populations
  4. Bounds: Partial identification when full identification isn’t possible
  5. Counterfactual dynamics bridge Dynamical and Observable: Use dynamic models to reason about alternative observable outcomes

21.8 Further Reading