---
title: "Getting started: state-space models with kalmix"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Getting started: state-space models with kalmix}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<!-- Render time: ~1 s under rmarkdown::render(); macOS arm64 (Darwin
     aarch64), R 4.5.2, one core. -->

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7,
                      fig.height = 4.5, dpi = 150)
set.seed(1)
has_ggplot2 <- requireNamespace("ggplot2", quietly = TRUE)
```

```{r library}
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()`](../reference/ssm.html) builds a linear-Gaussian state-space model
from its four system matrices and a prior on the initial state.
[`kalman_filter()`](../reference/kalman_filter.html) 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()`](../reference/rts_smoother.html) 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:

```{r local-level}
model <- ssm(
  transition = 1, observation = 1,
  state_cov = 0.01, obs_cov = 1,
  init_state = 0, init_cov = 10
)
model
```

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:

```{r filter}
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:

```{r smooth}
smoothed <- rts_smoother(fit)
rmse_smoothed <- rmse(smoothed@smoothed_mean[, 1], moisture)
```

```{r rmse-kable, echo = FALSE}
knitr::kable(
  data.frame(
    pass = c("filtered", "smoothed"),
    rmse = round(c(rmse_filtered, rmse_smoothed), 3)
  ),
  col.names = c("Pass", "RMSE against the true moisture index"),
  caption = "Root-mean-square error of each pass against the simulated truth."
)
```

```{r filter-plot, echo = has_ggplot2, eval = has_ggplot2, fig.cap = "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."}
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()
```

```{r filter-plot-skip, echo = FALSE, eval = !has_ggplot2, results = "asis"}
cat("ggplot2 is not installed on this build, so the filtered-and-smoothed",
    "figure is skipped.\n")
```

### 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:

```{r velocity}
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]
```

```{r velocity-kable, echo = FALSE}
knitr::kable(
  data.frame(
    quantity = "daily water-use rate (mm/day, synthetic)",
    truth = true_rate,
    recovered = round(recovered_rate, 3)
  ),
  col.names = c("Quantity", "True rate", "Recovered rate at day 150"),
  caption = paste(
    "The unobserved rate recovered purely from the noisy cumulative",
    "series."
  )
)
```

## Read

The smoothed pass has the lower RMSE of the two (`r round(rmse_smoothed, 3)`
against `r round(rmse_filtered, 3)` 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
`r round(recovered_rate, 3)` mm/day at the end of the series against a true
rate of `r true_rate` 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.

## What to read next

*Change-points and hidden states with kalmix* moves from a continuously
evolving state to abrupt shifts and discrete regimes. *Pairs trading with
kalmix* applies this same filter-and-smoother machinery to a mean-reverting
spread. *Assimilative causal inference with kalmix* reads a causal
information rate off the filter and smoother built here, grounded through a
model-adequacy check that can abstain.

## Reproduce

Seed `1` for both simulations.

```{r session-info}
sessionInfo()
```