# ADR-012: od-kalman — Davis & Nihan (1993) linear-Gaussian OD estimation from a time series of link counts

**Status:** accepted (implemented)
**Date:** 2026-07-08
**Deciders:** T2 estimation track — adding a time-series / covariance-aware baseline
**File:** `docs/design/adr-012-dn-kalman.md`

## Context

The six shipped T2 estimators (`prior`, `vzw-entropy`, `gls`, `spiess`, `spsa`,
`od-congested`) all collapse the observed link-count series to a single **mean**
vector (`counts.mean(axis=0)`) and fit that. They discard the two things Davis &
Nihan (1993, *Operations Research* 41(1):169-178, `davis1993large`) show carry
demand information: the **cross-link covariance** of the counts and their
**day-to-day autocorrelation**.

Davis & Nihan prove that as the traveler population grows, a general stochastic
assignment process converges to a deterministic mean (the SUE/UE loading, Prop 2)
plus a stationary linear-Gaussian VAR(1)/VARMA(1,1) fluctuation whose covariance
is their eq. (8) (Prop 3). Their stated application (§4) is exactly an estimator:
the VAR process *"admits a state-space description, which can then be combined
with the discrete-time Kalman filter to produce an algorithm for computing the
likelihood function of a time-series of observed link volume counts [...] permits
estimation of the underlying travel demands via prediction-error minimization."*

**A DN estimator cannot be built on the existing observation level.** `LinkCounts`
emits IID *per-link* Poisson counts (diagonal covariance) or exact repeats — no
cross-link and no temporal correlation. A DN estimator whitening a diagonal
covariance **is `gls`** (a rename the repo forbids; cf. the `od-congested`
"not-a-gls-rename" tests). The DN covariance only bites when the counts carry the
multinomial route structure and day-to-day persistence.

## Decision

1. **New observation level `DayToDayCounts`** (`observe/levels.py`, ADR-012)
   generating a benchmark realization of the DN Gaussian limit as a stationary
   VAR(1):
   `x(t) = x_UE + e(t)`, `e(t) = rho·e(t-1) + a(t)`, `a(t) ~ N(0, (1-rho^2) Q)`,
   so `e(t)` has stationary covariance `Q` and the series is **centered on the UE
   loading** `x_UE = P g`. `Q` is the route-exact DN multinomial covariance
   `Q = sum_j (g_j^2/N_j) Delta_j^T (diag p_j - p_j p_j^T) Delta_j`
   (`observe/_dn_process.py`), extracted from the equilibrium MSA/AON route sets
   so its link marginal equals the proportion column and the mean is exactly
   `x_UE`. `N_j = max(1, round(population_scale·g_j))`; the fluctuation vanishes as
   `1/population_scale` (Prop 2 SLLN) and reduces to `noise='none'` in the limit.
   `rho in [0,1)` is the day-to-day persistence dial; `rho=0` gives IID
   (cross-link-correlated only) counts.

2. **UE-centering (the load-bearing modeling choice).** The DN VARMA temporal
   structure is a *logit/SUE* phenomenon (the smooth `q_j(g)` feedback), but the
   T2 certifier pins **deterministic UE (bfw)** and the runner rejects SUE
   instances. Centering the VAR(1) on the UE loading with the DN innovation
   covariance keeps the **pinned-UE certificate valid with no SUE-certifier
   infra**, while preserving the two DN signals the estimator needs (spatial `Q`
   off-diagonals + temporal `rho`). It is honestly a *linearized Model-B process
   around UE* (Prop 3 form), not a re-simulation of the full logit adjustment.

3. **Estimator `od-kalman`** (`estimation/dn_kalman.py`,
   `DavisNihanKalmanEstimator`, paradigm `estimation`). For a **static** OD
   observed through a stationary Gaussian process the prediction-error / Kalman
   estimate is the GLS whitened by the covariance of the count **time-mean** — the
   steady-state Kalman gain for a constant state reduces to this batch BLUE. Two
   non-classical ingredients (`dn_gls_solve`):
   - **spatial:** whiten by the full (non-diagonal) sample covariance of the
     series, not `gls`'s diagonal Poisson `V`;
   - **temporal:** inflate the count-mean covariance by the AR(1) integrated
     autocorrelation `tau = (1+rho_hat)/(1-rho_hat)` (`ar1_tau`) — the
     effective-sample-size correction; `T` correlated observations carry the
     information of `T/tau` independent ones. `tau=1` recovers the naive mean.
   The solve runs the same congested outer fixed point as `gls`/`od-congested`
   (Cascetta & Postorino 2001) with the best-self-obs-RMSE safeguard (ADR-002
   Decision 3).

4. **Certificate (P1) unchanged.** The harness recomputes `od_feasible`,
   `od_rmse`, `obs_count_rmse` from the emitted OD matrix through the pinned bfw
   assignment (`metrics.estimation`); `self_report` is provenance only. No new
   scored key, no scenario-field or `FlowState`/`Trace`/`Evaluator` change — the
   golden Braess content hash `cf00f411…` is byte-identical (re-asserted in
   `tests/test_dn_kalman.py`).

5. **Analytic anchors** (`tests/test_dn_kalman.py`, all recomputed):
   - two-route DN covariance closed form `Var(link 0) = (D^2/N) p_A p_B`, with
     same-route `+` and cross-route `−` correlation (exact on link-disjoint
     routes);
   - single-sensor `dn_gls` closed form
     `g* = (g_pr/w^2 + p·c/s^2)/(1/w^2 + p^2/s^2)`;
   - `tau` recovers `(1+rho)/(1-rho)` on a DN series and `=1` when `rho=0`;
   - recovery of the planted `D` under full sensors at the finite-population floor.

## Alternatives considered

- **Ship on the existing IID Poisson counts (no new observation level):**
  rejected — the DN covariance is then diagonal and the estimator degenerates to
  `gls`; it would be a rename.
- **Full VARMA(1,1) marginal + a real logit/SUE adjustment process:** rejected for
  now — its mean is the SUE, which needs an SUE certifier (the runner rejects SUE
  T2 instances). The VAR(1)-around-UE slice is the minimal faithful realization
  that keeps the pinned-UE certificate; it also makes the pure-AR(1) form correct
  (no MA term to identify), sidestepping the VARMA MA-cancellation subtlety
  (p.175) — the `tau=(1+rho)/(1-rho)` constant is the *asymptotic* integrated
  autocorrelation, applied as `sigma_stat·tau/T` (a conservative finite-`T`
  effective-sample-size approximation). A full VARMA(1,1) + SUE-certified variant
  is a separately-scoped follow-up.
- **A sequential Kalman-filter recursion over the periods:** the batch
  DN-whitened GLS *is* the steady-state Kalman/BLUE for a static OD state and
  returns the same estimate; a sequential filter is more code for an identical
  result, so the established outer-loop/GLS machinery is reused.

## Consequences

The benchmark gains its **first time-series / covariance-aware** OD estimator: one
that exploits the cross-link and day-to-day structure of the counts rather than
their mean alone, with a sound, harness-recomputed certificate and hand-derived
anchors. All changes are additive; the golden Braess hash is provably preserved;
no output contract changed. Follow-ups: the full VARMA(1,1) + SUE-certified
variant, and reporting the DN posterior covariance of the estimate (not just the
point OD).

## Sourcing

Davis & Nihan (1993, *Operations Research* 41(1):169-178, `davis1993large`, JSTOR
stable/171951) is the primary and was **read** (the eq. (8) covariance recursion,
the VAR(1) logit form on p.175, the Prop 2 SUE mean, and the §4 Kalman /
prediction-error program). The route-level innovation covariance
`sum_j n_j [diag q_j - q_j q_j^T]` is COV[x|s] on p.171, written route-exactly;
it was cross-checked against an independent open implementation of the same eq.
(8) covariance. No number from the paper is reproduced — every anchor is a
hand-derived closed form.