Filtering and smoothing#

The Kalman filter#

The filter (DK §4.3) computes, for t = 1, …, n, the one-step prediction error and its variance,

\[v_t = y_t - d_t - Z_t a_t, \qquad F_t = Z_t P_t Z_t' + H_t,\]

the filtered state (DK eq. 4.24)

\[a_{t|t} = a_t + P_t Z_t' F_t^{-1} v_t, \qquad P_{t|t} = P_t - P_t Z_t' F_t^{-1} Z_t P_t,\]

and the prediction \(a_{t+1} = c_t + T_t a_{t|t}\), \(P_{t+1} = T_t P_{t|t} T_t' + R_t Q_t R_t'\). The log likelihood is (DK eq. 7.2)

\[\log L = -\frac{np}{2} \log 2\pi - \frac12 \sum_t \big(\log|F_t| + v_t' F_t^{-1} v_t\big).\]

The examples on this page run on a Representation with known matrices. For a model, mod.filter(params) and mod.smooth(params) do the same at given parameters, res.filter() and res.smooth() at the estimates, and mod.representation(params) returns the representation for the other methods below.

>>> y = np.loadtxt("data/nile.csv", delimiter=",", skiprows=1)[:, 1]
>>> rep = ss.Representation(y, k_states=1)
>>> rep["design"] = rep["transition"] = rep["selection"] = [[1.0]]
>>> rep["obs_cov"], rep["state_cov"] = [[15099.0]], [[1469.1]]
>>> rep.initialize_diffuse()
>>> f = rep.filter()
>>> f.nobs_diffuse, round(f.llf, 4)
(1, -633.4646)

loglike() returns the same number without storing the filter output; estimation uses it.

Univariate treatment#

With filter_method = FILTER_UNIVARIATE the filter processes the elements of \(y_t\) one at a time (DK §6.4). It needs no matrix inverse and handles singular \(F_t\). When \(H_t\) is not diagonal, each period’s observed block is first factorized, \(H_{oo} = L D L'\), and the observation equation transformed by \(L^{-1}\) (DK §6.4.3); this works with time-varying H and any pattern of missing values.

The reported \(v_t\) and \(F_t\) are in the original coordinates in every period; see Output coordinates of the univariate treatment.

Diffuse periods#

With diffuse states the exact initial filter (DK §5.2) runs until \(P_\infty\) vanishes, after nobs_diffuse periods. By default it processes those periods element by element (DK §5.2.5). diffuse_method = DIFFUSE_MULTIVARIATE uses the multivariate form where \(F_\infty\) is zero or nonsingular. The per-period choice is in the filter output’s method array (Fortran).

Steady state#

In a time-invariant model \(P_t\) converges. Once the relative change falls below tol_steady (default 1e-15), the filter holds \(P_t\), \(F_t\) and \(K_t\) fixed and updates only the state means (DK §4.3.4), until a missing observation. t_steady reports where this began. steady_state() computes the limit directly:

>>> P, F = rep.steady_state()
>>> round(float(F[0, 0]), 2)
20600.26

See The steady state.

State and disturbance smoothing#

The smoother (DK §4.4-4.5) runs backwards through

\[r_{t-1} = Z_t' F_t^{-1} v_t + L_t' r_t, \qquad N_{t-1} = Z_t' F_t^{-1} Z_t + L_t' N_t L_t,\]

with \(L_t = T_t - K_t Z_t\), and gives the smoothed state \(\hat\alpha_t = a_t + P_t r_{t-1}\) (DK eq. 4.39), its variance \(V_t = P_t - P_t N_{t-1} P_t\) (DK eq. 4.43), and the smoothed disturbances with their variances.

>>> sm = rep.smooth()
>>> sm.smoothed_state.shape, sm.smoothed_state_cov.shape
((1, 100), (1, 1, 100))

Other smoothers#

DK §4.4-4.8 give further algorithms; each is a method of Representation:

fast_smoother() (DK §4.6.2): means only;
classical_smoother() (DK §4.6.1), two_filter_smoother() (DK §4.6.4), whittle_smoother() (DK §4.6.3);
fixed_point_smoother() and fixed_lag_smoother() (DK §4.4.6), and update_smoothed() (DK §4.4.5) for observations that arrive later;
smoothed_state_cov_between() and smoothed_state_autocov() (DK §4.7);
filtered_state_weights() and smoothed_state_weights() (DK §4.8);
smooth(method="sqrt"), the square root filter and smoother (DK §6.3).

They agree with the standard smoother, which the tests check; the Whittle recursion loses accuracy on long series, as DK note.

>>> np.allclose(rep.fast_smoother(), sm.smoothed_state)
True