Assimilative causal inference with kalmix

Causality as an assimilation problem

The filter estimates the state from the past alone, giving \(N(x_{t\mid t}, P_{t\mid t})\); the smoother adds the future, giving \(N(x_{t\mid n}, P_{t\mid n})\). The causal information the measurement carries about the state at \(t\) is the relative entropy (Kullback-Leibler divergence) of the smoother from the filter,

\[I_t = D_{\mathrm{KL}}\!\left( N(x_{t\mid n}, P_{t\mid n}) \,\Vert\, N(x_{t\mid t}, P_{t\mid t}) \right).\]

For two Gaussians this is a closed form, which kalmix computes natively. Time averaged, and integrated over an expanding future window, the divergence yields a threshold-free causal information rate, a decision lead-time in the time units of the series.

A known-truth test: a recharge oscillator

The cleanest way to trust a causal read-out is to run it on a system whose truth is controlled. We simulate a synthetic series from a noise-sustained recharge oscillator (Jin, 1997), the canonical conceptual model of the El Nino-Southern Oscillation (ENSO). The hidden state \((T, h)\) is a sea-surface-temperature anomaly \(T\) and a thermocline-depth anomaly \(h\); their coupled, weakly damped dynamics produce an irregular oscillation of roughly four years. We observe only a scalar that warm \(T\) drives through a coupling kappa, standing in for a crop-water-stress index that rises in an El Nino year. When kappa is zero the observation tells us nothing about the hidden cause.

simulate_enso <- function(n = 600, kappa = 0.6, dt = 1 / 12, obs_sd = 0.5) {
  drift <- matrix(c(0, 1.6, -1.6, -0.2), nrow = 2, byrow = TRUE)
  amat <- diag(2) + drift * dt
  z <- matrix(0, nrow = n, ncol = 2)
  x <- numeric(n)
  for (t in 2:n) {
    z[t, ] <- as.numeric(amat %*% z[t - 1, ]) +
      c(rnorm(1, sd = 0.35 * sqrt(dt)), rnorm(1, sd = 0.25 * sqrt(dt)))
    x[t] <- kappa * z[t, 1] + rnorm(1, sd = obs_sd)
  }
  list(x = x, z = z, dt = dt, amat = amat)
}

sim <- simulate_enso(kappa = 0.6)

We give aci() a state-space model matching the data-generating process. The observation operator maps the first hidden component (the SST anomaly) to the observed series through the coupling:

model <- ssm(
  transition = sim$amat,
  observation = matrix(c(0.6, 0), nrow = 1),
  state_cov = diag(c(0.35^2 * sim$dt, 0.25^2 * sim$dt)),
  obs_cov = 0.5^2,
  init_state = c(0, 0),
  init_cov = diag(c(1, 1))
)

fit <- aci(model, sim$x, dt = sim$dt)
fit
#> <aci_fit>: assimilative causal inference (600 steps)
#>   mean causal information : 0.6204 nats
#>   decision lead-time      : 0.527 time units
#>   grounding               : ungrounded (mechanism_unverified)

recovery_cor <- cor(fit@smoother@smoothed_mean[, 1], sim$z[, 1])

Two things to check against the known truth. First, the smoother recovers the hidden SST anomaly from the observed series alone, the inverse-problem essence, at a correlation of 0.983 with the simulated truth. Second, the causal information series should be positive where the future of the series is informative about the hidden state.

years <- (seq_along(sim$x) - 1) * sim$dt
plot_data <- rbind(
  data.frame(year = years, value = sim$z[, 1], series = "truth (simulated)",
             panel = "Hidden ENSO state"),
  data.frame(year = years, value = fit@smoother@smoothed_mean[, 1],
             series = "smoother estimate", panel = "Hidden ENSO state"),
  data.frame(year = years, value = fit@causal_information,
             series = "causal information",
             panel = "Causal information I_t (nats)")
)
panel_levels <- c("Hidden ENSO state", "Causal information I_t (nats)")
plot_data$panel <- factor(plot_data$panel, levels = panel_levels)

ggplot2::ggplot(plot_data, ggplot2::aes(x = year, y = value, colour = series)) +
  ggplot2::geom_line(linewidth = 0.5) +
  ggplot2::facet_wrap(~panel, ncol = 1, scales = "free_y") +
  ggplot2::scale_colour_viridis_d(name = NULL, end = 0.8) +
  ggplot2::labs(
    x = "year", y = NULL,
    title = "Recovered state and causal information, coupled system"
  ) +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "bottom")

Top: the hidden ENSO state recovered from the observed index alone, against the truth used to simulate it. Bottom: the per-step causal information series that recovery is built from.

The smoother trace tracks the truth closely in the top panel although the model never observes the state directly, and the causal-information trace in the bottom panel rises and falls with the oscillation rather than sitting flat at zero, which is the discriminating signature the null case below is checked against.

The null carries no causal information

Under the null of no coupling between the hidden state and the observation, the future of the series adds nothing the filter did not already have, so the causal information collapses to zero and the lead-time with it:

sim_null <- simulate_enso(kappa = 0)
model_null <- ssm(
  transition = sim_null$amat,
  observation = matrix(c(1e-6, 0), nrow = 1),
  state_cov = diag(c(0.35^2 * sim_null$dt, 0.25^2 * sim_null$dt)),
  obs_cov = 0.5^2,
  init_state = c(0, 0),
  init_cov = diag(c(1, 1))
)
fit_null <- aci(model_null, sim_null$x, dt = sim_null$dt)

A divergence between two Gaussians, cross-checked

The metric aci() reads off the two passes is the relative entropy between the smoother and filter Gaussians at each step. This is canonical Gaussian arithmetic, so kalmix implements it natively and stays self-contained. Where the kernR package is installed, its relative_entropy() provides an independent check on the value: the package’s test suite asserts the two agree and skips the comparison when kernR is absent, so kalmix builds and checks cleanly either way.

Grounding the verdict: abstention under model error

ACI reads its causal verdict off a state-space model, so the verdict is only as good as the model. A model that cannot even describe the observed series cannot support a causal claim about it. aci() therefore grounds every read-out through innovation_diagnostics(), run on the standardised one-step-ahead innovations. When the innovations reject the model it is declared inadequate and the verdict abstains, with a distinct grounding token instead of a confident lead-time:

## A frozen-state local level cannot track the oscillating ENSO index: its
## innovations are autocorrelated and mis-scaled, so ACI abstains.
frozen <- ssm(transition = 1, observation = 1, state_cov = 1e-6, obs_cov = 1,
              init_state = 0, init_cov = 10)
abstained <- aci(frozen, sim$x, dt = sim$dt)
abstained@grounding
#> [1] "ungrounded"
abstained@grounding_reason
#> [1] "model_inadequate"
abstained@adequacy@adequate
#> [1] FALSE
abstention_classes <- class(abstained)
abstention_classes
#> [1] "orchestra_refusal" "kalmix_abstention" "kalmix::aci_fit"  
#> [4] "S7_object"

## Adequate but unverified by default; grounded only with dated provenance.
adequate_unverified <- aci(model, sim$x, dt = sim$dt)@grounding
adequate_grounded <- aci(
  model, sim$x, dt = sim$dt,
  mechanism = list(source = "validated against the OSSE truth",
                    verified_on = as.Date("2026-06-14"))
)@grounding

The lead-time without re-smoothing

The decision lead-time integrates a divergence profile over an expanding future window. Formed directly, each window length reruns the smoother on a truncated series. kalmix instead grows the smoothed estimate of the anchor state forward in place as each later observation is assimilated (a fixed-point smoother recursion that accumulates the Rauch-Tung-Striebel gains), so the whole profile costs one filter-smoother pass and a light forward recursion per anchor, with no re-smoothing. This online engine is the default; the expanding engine is retained as an independent cross-check:

cir_online <- causal_information_rate(model, sim$x, dt = sim$dt,
                                       engine = "online")
cir_expanding <- causal_information_rate(model, sim$x, dt = sim$dt,
                                          engine = "expanding")

Because the online recursion never truncates the series, it also applies to time-varying state-space models, where the expanding construction cannot: a model whose system matrices change step by step still yields a lead-time from a single forward pass. The two engines agreeing shows only that the smoother recursion is implemented consistently between them: they share the quadrature and the reduction, so their agreement is not evidence about the range functional itself. That functional is graded separately, in the package’s own conformance suite (tests/testthat/test-cir-conformance.R), which reduces acir‘s own divergence profile with kalmix’s reducer and requires the objective and subjective ranges of Andreou, Chen and Bollt (2026) to come back; acir is in turn graded cell-by-cell against the method authors’ MATLAB reference. That chain, back to the authors’ own code, is what makes it an oracle here, and it is not recomputed in this vignette.

A real ENSO index

The known-truth oscillator proves the read-out is internally consistent; the package’s test suite also runs aci() on the real NOAA Oceanic Nino Index (ONI) as an external check, in tests/testthat/test-aci.R, on a fixture shipped with the package (tests/testthat/fixtures/oni_noaa.csv, the CPC ONI series, 1980 to 2020, public-domain US Government data, NOAA Climate Prediction Center). The same model and data are reproduced here directly, rather than cited unrun, so the number below is computed in this vignette:

oni_path <- file.path("..", "tests", "testthat", "fixtures", "oni_noaa.csv")
has_oni <- file.exists(oni_path)

oni <- read.csv(oni_path, stringsAsFactors = FALSE)$oni

## A two-state local oscillator: a damped rotation captures the quasi-
## periodic ENSO swing, observed through the index with measurement noise.
dt_year <- 1 / 12
angle <- 2 * pi * dt_year / 3.8
rho <- 0.96
rot <- rho * matrix(
  c(cos(angle), -sin(angle), sin(angle), cos(angle)),
  nrow = 2L, byrow = TRUE
)
oni_model <- ssm(
  transition = rot,
  observation = matrix(c(1, 0), nrow = 1L),
  state_cov = diag(c(0.05, 0.05)),
  obs_cov = 0.2,
  init_state = c(oni[1L], 0),
  init_cov = diag(c(1, 1))
)
oni_fit <- aci(oni_model, oni, dt = dt_year)
oni_tracked <- cor(oni_fit@smoother@smoothed_mean[, 1L], oni)

Two things bound how far that figure should be pushed. The quadrature is evaluated at every lag in the window rather than at a dozen points, and warns when a grid the caller thins has not converged, so the number is not grid-dependent at the default. And an anchor whose future window runs past the end of the record is right-censored: its range is integrated over the lags the record supports and returned as a lower bound with censored set, rather than extrapolated. A rolling read-out near the end of a record is censored by construction, and a censored lead-time cannot carry a grounded verdict.