Getting started: state-space models with kalmix

library(kalmix)

Why

A soil-moisture probe in one paddock logs a noisy reading every day of the season. What is the moisture level actually doing underneath the sensor noise, and could a quantity that was never measured directly, such as how fast the level is changing, be read off that noisy trace alone? A state-space model answers exactly this: a latent state evolves on its own dynamics and is seen only through noisy observations, and the recursions in kalmix recover it. This vignette works the recipe on two synthetic soil-moisture and crop-water-use series, small enough to see the whole mechanism at once. The Change-points and hidden states with kalmix vignette covers discrete-state models, and Pairs trading with kalmix covers the mean-reverting-spread application built on the same machinery.

What

ssm() builds a linear-Gaussian state-space model from its four system matrices and a prior on the initial state. kalman_filter() runs the forward recursion of Kalman (1960): at each step it predicts the state forward and then corrects it with the new observation, conditioning only on the past. rts_smoother() runs the backward pass of Rauch, Tung and Striebel (1965) over a completed filter, conditioning every state estimate on the whole series instead of the past alone.

Do

A linear-Gaussian state-space model

The general model has a latent state \(x_t\) that evolves through a transition matrix and Gaussian process noise, observed through an observation matrix and Gaussian observation noise:

\[x_t = A_t\, x_{t-1} + w_t, \qquad w_t \sim N(0, Q_t),\] \[y_t = B_t\, x_t + v_t, \qquad v_t \sim N(0, R_t).\]

The simplest case is the local level (a random walk seen through observation noise), which can be written with bare scalars. We take it to represent a soil-moisture index that drifts slowly over a season, in arbitrary synthetic units:

model <- ssm(
  transition = 1, observation = 1,
  state_cov = 0.01, obs_cov = 1,
  init_state = 0, init_cov = 10
)
model
#> <ssm>: linear-Gaussian state-space model (state 1, obs 1)
#>   transition  : static
#>   observation : static
#>   state_cov   : static
#>   obs_cov     : static

The print method above shows the model exactly as built: a one-dimensional state, unit transition and observation, and the process and observation noise variances that make the observation channel ten times noisier than the state’s own drift.

Filtering: recovering the state online

We simulate a slowly-drifting soil-moisture index seen through instrument noise, then recover it with the forward recursion:

moisture <- cumsum(rnorm(150, sd = 0.1))
y <- moisture + rnorm(150)

fit <- kalman_filter(model, y)
rmse <- function(a, b) sqrt(mean((a - b)^2))
rmse_filtered <- rmse(fit@filtered_mean[, 1], moisture)

Smoothing: using the whole series

The filter conditions only on the past; the smoother adds the future, which is strictly more informative:

smoothed <- rts_smoother(fit)
rmse_smoothed <- rmse(smoothed@smoothed_mean[, 1], moisture)
Root-mean-square error of each pass against the simulated truth.
Pass RMSE against the true moisture index
filtered 0.269
smoothed 0.198 
band_data <- data.frame(
  day = rep(seq_len(150), 2),
  estimate = c(fit@filtered_mean[, 1], smoothed@smoothed_mean[, 1]),
  se = c(
    sqrt(vapply(fit@filtered_cov, function(p) p[1, 1], numeric(1))),
    sqrt(vapply(smoothed@smoothed_cov, function(p) p[1, 1], numeric(1)))
  ),
  pass = rep(c("filtered", "smoothed"), each = 150)
)
band_data$lower <- band_data$estimate - 1.96 * band_data$se
band_data$upper <- band_data$estimate + 1.96 * band_data$se

obs_data <- data.frame(day = seq_len(150), observed = y, truth = moisture)

ggplot2::ggplot() +
  ggplot2::geom_point(
    data = obs_data, ggplot2::aes(x = day, y = observed),
    colour = "grey60", size = 0.8
  ) +
  ggplot2::geom_line(
    data = obs_data, ggplot2::aes(x = day, y = truth),
    linetype = "dashed", colour = "black", linewidth = 0.4
  ) +
  ggplot2::geom_ribbon(
    data = band_data,
    ggplot2::aes(x = day, ymin = lower, ymax = upper, fill = pass),
    alpha = 0.22
  ) +
  ggplot2::geom_line(
    data = band_data, ggplot2::aes(x = day, y = estimate, colour = pass),
    linewidth = 0.7
  ) +
  ggplot2::scale_colour_viridis_d(name = "pass", end = 0.8) +
  ggplot2::scale_fill_viridis_d(name = "pass", end = 0.8) +
  ggplot2::labs(
    x = "day of season", y = "soil-moisture index (synthetic units)",
    title = "Filtered and smoothed soil-moisture estimate"
  ) +
  ggplot2::theme_minimal()
The noisy daily readings (grey points) against the true soil-moisture index (dashed line), with the filtered and smoothed estimates and their 95 per cent bands.

The noisy daily readings (grey points) against the true soil-moisture index (dashed line), with the filtered and smoothed estimates and their 95 per cent bands.

A vector state

The state can have more than one dimension. A constant-rate model tracks a cumulative quantity and an unobserved rate; here, cumulative crop water use and the unobserved daily water-use rate that drives it. Only the cumulative total is measured, and the model recovers the hidden rate from how that total changes:

velocity_model <- ssm(
  transition = matrix(c(1, 0, 1, 1), nrow = 2),
  observation = matrix(c(1, 0), nrow = 1),
  state_cov = diag(c(1e-5, 1e-5)),
  obs_cov = 1,
  init_state = c(0, 1),
  init_cov = diag(c(10, 10))
)

true_rate <- 0.5
cum_use <- true_rate * seq_len(150) + rnorm(150)
vfit <- kalman_filter(velocity_model, cum_use)
recovered_rate <- vfit@filtered_mean[150, 2]
The unobserved rate recovered purely from the noisy cumulative series.
Quantity True rate Recovered rate at day 150
daily water-use rate (mm/day, synthetic) 0.5 0.485

Read

The smoothed pass has the lower RMSE of the two (0.198 against 0.269 for the filter), as it must: it has seen the whole series rather than the past alone. The figure shows why the two differ visually as well as numerically: the blue-green filtered band sits close to the observations early in the record, before much of the series has arrived to correct it, while the smoothed band is narrower and tighter to the dashed truth line throughout, because every point on it has already seen the observations that come after it.

In the vector-state example, the daily water-use rate was never observed; only the cumulative total was. The filter still recovers a rate of 0.485 mm/day at the end of the series against a true rate of 0.5 mm/day, inferred entirely from the curvature of the noisy cumulative trace.

Limits

Both examples here are synthetic, chosen to keep the recursions fast and the truth known, not to represent a validated soil-moisture or water-use record. The system matrices are supplied by the caller and taken as given: this vignette does not check whether a model is an adequate description of the series it is applied to, nor does it estimate unknown noise levels from data, which ssm_fit() does (see The Elliott-van der Hoek-Malcolm spread model with kalmix). Every model shown here is linear-Gaussian and single-regime; a series that switches between discrete states or regimes needs the machinery in the next vignette. Grounding a causal claim against a model’s own adequacy diagnostics, including the abstention a diagnostic failure produces, is shown in Assimilative causal inference with kalmix, not here.

Reproduce

Seed 1 for both simulations.

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] kalmix_0.6.0.9000 rmarkdown_2.32   
#> 
#> loaded via a namespace (and not attached):
#>  [1] vctrs_0.7.3        cli_3.6.6          knitr_1.52         rlang_1.3.0       
#>  [5] xfun_0.61          otel_0.2.0         S7_0.2.2           jsonlite_2.0.0    
#>  [9] labeling_0.4.3     glue_1.8.1         buildtools_1.0.0   htmltools_0.5.9   
#> [13] maketools_1.3.2    sys_3.4.3          sass_0.4.10        scales_1.4.0      
#> [17] grid_4.6.1         evaluate_1.0.5     jquerylib_0.1.4    fastmap_1.2.0     
#> [21] yaml_2.3.12        lifecycle_1.0.5    compiler_4.6.1     RColorBrewer_1.1-3
#> [25] farver_2.1.2       digest_0.6.39      viridisLite_0.4.3  R6_2.6.1          
#> [29] bslib_0.12.0       withr_3.0.3        tools_4.6.1        gtable_0.3.6      
#> [33] ggplot2_4.0.3      cachem_1.1.0