Model specification (same as dlmGenSys)
Observation times (length n). Timesteps Δt_k = t_{k} - t_{k-1} are computed from consecutive differences. The first Δt is t[1]-t[0] (t[0] is treated as the baseline; the first transition is from prior to t[0] so Δt[0] is set to 1.0, matching the uniform-timestep convention).
Process noise std devs (diagonal of √Q_c per unit time). Length determines which states have noise (same as DlmFitOptions.processStd).
Time-varying system matrices { G[n,m,m], W[n,m,m], F[m], m }
Generate time-varying DLM system matrices G(Δt) and W(Δt) for non-uniform timesteps in closed form.
G(Δt) is the genuine continuous-time discretization: G_k = expm(F_c · Δt_k) where F_c is the continuous-time state matrix. For polynomial trend this evaluates to the Jordan block with entries Δt^j/j! on the j-th superdiagonal; no numerical matrix exponential is needed.
W(Δt) is the discrete-time noise accumulation sum, analytically continued to non-integer Δt via Faulhaber sums: W_k = Σ_{j=0}^{Δt-1} G^j · W_1 · (G^j)ᵀ Because G is a Jordan block, G^j has polynomial entries in j and the sum reduces to closed-form Faulhaber sums. Evaluating those polynomials at non-integer Δt is an analytic continuation of the discrete-time formula. This is NOT the continuous-time SDE integral ∫ expm(F_c·s) Q_c expm(F_c·s)' ds. The Faulhaber formula is the correct choice when Δt represents "skipping Δt discrete time steps" rather than "evolving a continuous SDE for Δt time units". It collapses exactly to W_1 = diag(w²) at unit spacing.
Supported components: polynomial trend (order 0, 1, 2), trigonometric harmonics, and AR components (integer timestep spacing only). Throws if fullSeasonal is requested (purely discrete-time construct with no natural continuous-time extension). AR components require all Δt to be positive integers; G_AR^d is computed via binary matrix exponentiation and W_AR(d) = Σ_{k=0}^{d-1} C^k · W₁ · (C^k)ᵀ by direct summation.