ssfortran.Representation#

class ssfortran.Representation(endog, k_states, k_posdef=None)#

A linear Gaussian state space model and its algorithms.

The system matrices start at zero, except selection, which starts as the leading m x r identity; they are set by item assignment. See the module description for their names and shapes. The initial state must be set with one of the initialize_* methods before filtering.

Parameters:

endogarray_like, shape (n,) or (n, p): Observations, time first as in statsmodels. NaN marks a missing value; any subset of a period’s elements may be missing.
k_statesint: Number of states m.
k_posdefint, optional: Number of state disturbances r. The default is k_states.

Attributes:

nobs, k_endog, k_states, k_posdefint: n, p, m and r.
filter_methodint: FILTER_CONVENTIONAL (default) or FILTER_UNIVARIATE. Set-only.
diffuse_methodint: DIFFUSE_UNIVARIATE (default) or DIFFUSE_MULTIVARIATE. Set-only.
loglikelihood_burnint: Number of leading periods left out of the log likelihood. Set-only.
marginal_likelihoodbool: Report the marginal instead of the diffuse log likelihood (DK §7.2.6). Set-only.
tol_diffusefloat: Threshold on \(F_\infty\) and \(\|P_\infty\|_F^2\) in the diffuse periods; default 1e-10, as in statsmodels. Set-only.
tol_steadyfloat: Relative change in \(P_t\) below which the filter holds P, F and K fixed (DK §4.3.4); default 1e-15. Negative turns the shortcut off. Set-only.

See also

StructuralModel, MappedModel, MLEModel: Models with parameters.

Notes

The filter treats diffuse periods element by element (DK §5.2.5, §6.4) and reports \(v_t\) and \(F_t\) in the original coordinates; the per-element quantities stay internal. With time-invariant system matrices the conventional filter stops updating \(P_t\) once it converges and resumes after a missing observation.

Examples

The local level model (DK ch. 2) with exact diffuse initialization:

>>> import numpy as np
>>> import ssfortran as ss
>>> rng = np.random.default_rng(0)
>>> y = rng.standard_normal(100).cumsum() + rng.standard_normal(100)
>>> rep = ss.Representation(y, k_states=1)
>>> rep["design"] = rep["transition"] = rep["selection"] = [[1.0]]
>>> rep["obs_cov"] = [[1.0]]
>>> rep["state_cov"] = [[1.0]]
>>> rep.initialize_diffuse()
>>> round(rep.loglike(), 4)
-189.1596
>>> res = rep.smooth()
>>> res.smoothed_state.shape
(1, 100)

__getitem__(name)#: Copy of a system matrix; time-invariant ones without the time axis.

__setitem__(name, value)#: Set a system matrix, or endog with shape (p, n).

add_state_restrictions(restriction, value)#

Impose linear restrictions on the states (DK §6.6).

The restrictions \(R^*_t \alpha_t = r^*_t\) enter as exact observations.

Parameters:

restrictionarray_like, shape (q, m) or (q, m, n): \(R^*_t\).
valuearray_like, shape (q, n): \(r^*_t\); NaN where a restriction does not apply.

Returns:

Representation: The model with the restrictions as extra observations. Filter it with filter_method = FILTER_UNIVARIATE: once a restricted direction has no state noise, \(F_t\) is singular.

augmented()#

Run the augmented Kalman filter and smoother (DK §5.7).

Returns:

AugmentedResults: The filter, the estimate of the diffuse elements and three log likelihoods.

See also

marginal_likelihood: Option for the marginal log likelihood.

auxiliary_residuals(vector=False)#

Compute the standardized smoothed disturbances (DK §7.5).

Large values point to outliers (observation disturbances) and structural breaks (state disturbances).

Parameters:

vectorbool, optional: Standardize each period’s disturbance vector by the Cholesky factor of its variance instead of element by element.

Returns:

epsndarray, shape (p, n): Auxiliary residuals of the observation equation.
etandarray, shape (r, n): Auxiliary residuals of the state equation. NaN where the variance is zero.

classical_smoother()#

Compute the smoothed state by the classical fixed-interval smoother.

The Rauch-Tung-Striebel form (DK §4.6.1), from the filtered states. Not available with diffuse states.

Returns:

alphahatndarray, shape (m, n): Smoothed state.
Vndarray, shape (m, m, n): Its variance.

collapse()#

Collapse the observations to the state dimension (DK §6.5).

Each \(y_t\) is replaced by its generalized least squares projection on \(\alpha_t\), which gives the same filtered and smoothed states with p reduced to m (Jungbacker and Koopman 2015).

Returns:

collapsedRepresentation: The collapsed model.
adjustndarray, shape (n,): Log likelihood adjustment: self.loglike() == collapsed.loglike() + adjust.sum().

Raises:

StateSpaceError: With code 2 if a period’s \(Z_t\) (observed rows) lacks full column rank.

de_jong_penzer()#

Compute the shock statistics of de Jong and Penzer (1998) (DK §7.5).

Returns:

statendarray, shape (m, n): \(r_{i,t} / \sqrt{N_{ii,t}}\).
observationndarray, shape (p, n): \(e_{i,t} / \sqrt{D_{ii,t}}\). NaN where undefined.

djs_simulation_smoother(rng=None, variates=None)#

Draw states and disturbances by the de Jong-Shephard method.

DK §4.9.3: the disturbances are drawn backwards from their distribution given the data and the later draws.

Parameters:

rngnumpy.random.Generator or int, optional: Source of the standard normal variates.
variatestuple of array_like, optional: (u_init (m,), u_eps (p, n), u_eta (r, n)); overrides rng.

Returns:

alphandarray, shape (m, n): Draw of the states.
epsndarray, shape (p, n): Draw of the observation disturbances.
etandarray, shape (r, n): Draw of the state disturbances.

See also

simulation_smoother: The mean-correction method (DK §4.9.2).

em(maxiter=500, tol=1e-08, diagonal_obs_cov=False, diagonal_state_cov=False)#

Estimate H and Q by the EM algorithm (DK §7.3.4).

Updates obs_cov and state_cov in place. The E-step uses the smoothed disturbances and their variances; for missing elements these are the conditional moments given the observed elements.

Parameters:

maxiterint, optional: Maximum number of iterations.
tolfloat, optional: Stop when the log likelihood improves by less than tol.
diagonal_obs_cov, diagonal_state_covbool, optional: Restrict H or Q to be diagonal.

Returns:

llffloat: Log likelihood at the start of the last iteration.
niterint: Number of iterations.
llf_pathndarray, shape (niter,): Log likelihood at each iteration; never decreasing.

Notes

H and Q must be time-invariant, the initialization must not depend on them (not stationary), and there must be no burn-in.

property endog#: Observations, shape (p, n).

fast_smoother()#

Compute the smoothed state by the fast state smoother (DK §4.6.2).

A backward pass for \(r_t\) only, then the forward recursion \(\hat\alpha_{t+1} = c_t + T_t \hat\alpha_t + R_t Q_t R_t' r_t\); no variances.

Returns:

ndarray, shape (m, n): Smoothed state, equal to that of smooth.

filter(method='conventional')#

Run the Kalman filter.

Parameters:

method{“conventional”, “sqrt”}, optional: "sqrt" propagates square roots of \(P_t\) (DK §6.3); it gives the same output and does not support diffuse states.

Returns:

FilterResults: The filter output.

See also

smooth: Filter and smoother.
loglike: Log likelihood without storing the output.

filtered_state_weights()#

Compute the weights of past observations in the filtered state.

DK §4.8.2; intercepts and the initial mean aside.

Returns:

Wandarray, shape (m, p, n, n): Wa[:, :, t, j] is the weight of \(y_j\) in \(a_t\).
Wattndarray, shape (m, p, n, n): The same for \(a_{t|t}\).

fixed_lag_smoother(lag)#

Estimate each state from the observations up to lag periods later.

DK §4.4.6.

Parameters:

lagint: Number of later observations used.

Returns:

alphahatndarray, shape (m, n): \(E(\alpha_s | y_1, \dots, y_{s+lag})\); NaN in the last lag periods.
Vndarray, shape (m, m, n): Its variance.

fixed_point_smoother(t)#

Track the estimate of one state as observations arrive (DK §4.4.6).

Parameters:

tint: The period (0-based).

Returns:

meanndarray, shape (m, n - t): Column k is \(E(\alpha_t | y_1, \dots, y_{t+k})\).
varndarray, shape (m, m, n - t): The corresponding variances.

forecast(steps)#

Forecast the observations past the end of the sample (DK §4.11).

Parameters:

stepsint: Number of periods ahead.

Returns:

meanndarray, shape (p, steps): \(E(y_{n+j} | Y_n)\).
covndarray, shape (p, p, steps): \(Var(y_{n+j} | Y_n)\).

Raises:

StateSpaceError: With code 4 if a system matrix varies over time. Append NaN observations, with the future matrices, to forecast such models.

initialize(a1, Pstar, Pinf)#

Initialize with \(\alpha_1 \sim N(a_1, P_* + \kappa P_\infty)\).

The general form of DK §5.1, with \(\kappa \to \infty\) handled exactly.

Parameters:

a1array_like, shape (m,): Mean.
Pstararray_like, shape (m, m): Variance of the proper part.
Pinfarray_like, shape (m, m): Variance direction of the diffuse part; usually a selection \(A A'\).

initialize_approximate_diffuse(variance=1000000.0)#

Initialize with mean zero and a large variance, variance * I.

Parameters:

variancefloat, optional: Variance of each state; default 1e6, as in statsmodels.

See also

initialize_diffuse: The exact treatment (DK §5.2).

initialize_block(start, stop, kind, a1=None, P1=None, variance=1000000.0)#

Initialize a block of states, leaving the others as they are.

Blocks combine, for example a diffuse trend with a stationary ARMA block (DK §5.6).

Parameters:

start, stopint: The block is states start:stop (0-based, like a slice).
kindint: INIT_KNOWN, INIT_APPROX_DIFFUSE, INIT_STATIONARY or INIT_DIFFUSE.
a1, P1array_like, optional: Mean (stop - start,) and variance for INIT_KNOWN.
variancefloat, optional: Variance for INIT_APPROX_DIFFUSE.

Notes

A stationary block must not depend on states outside it through T.

Examples

A diffuse level with a stationary AR(1) disturbance:

>>> rep = ss.Representation(np.zeros(20), k_states=2)
>>> rep["transition"] = [[1.0, 0.0], [0.0, 0.5]]
>>> rep.initialize_block(0, 1, ss.INIT_DIFFUSE)
>>> rep.initialize_block(1, 2, ss.INIT_STATIONARY)

initialize_diffuse()#

Initialize every state as exact diffuse (DK §5.2).

The initial variance is \(\kappa I\) with \(\kappa \to \infty\), handled exactly by the exact initial Kalman filter; the log likelihood is then the diffuse log likelihood (DK §7.2.2).

See also

initialize_approximate_diffuse: A large finite variance instead.

initialize_known(a1, P1)#

Initialize with a known mean and variance.

Parameters:

a1array_like, shape (m,): Mean of the initial state.
P1array_like, shape (m, m): Variance of the initial state.

See also

initialize_diffuse: Exact diffuse initialization.
initialize_stationary: The unconditional distribution.
initialize_block: Different initializations for blocks of states.

initialize_stationary()#

Initialize with the unconditional distribution of the state.

The mean is \((I - T)^{-1} c\) and the variance solves \(P = T P T' + R Q R'\) (DK §5.6.2), recomputed from the current matrices at every filter run. The model must be time-invariant in T, R, Q and c at the first period.

Raises:

StateSpaceError: From the filter, with code 6, if T has an eigenvalue on or outside the unit circle.

innovation_transition(t)#

Compute \(L_t\), the transition of the state prediction error.

\(a_{t+1} - \alpha_{t+1} = L_t (a_t - \alpha_t) + \dots\) with \(L_t = T_t - K_t Z_t\) (DK §4.3).

Parameters:

tint: Period (0-based), past the diffuse periods.

Returns:

ndarray, shape (m, m): \(L_t\).

least_squares_residuals(start, stop)#

Compute least squares residuals for regression effects (DK §6.2.4).

Parameters:

start, stopint: The regression coefficients are states start:stop (0-based); they must be diffuse and constant.

Returns:

ndarray, shape (p, n): Residuals with the coefficients at their full-sample estimate.

loglike()#

Evaluate the log likelihood without storing the filter output.

Returns:

float: The log likelihood (DK eq. 7.2), or the diffuse log likelihood (DK §7.2.2) with diffuse states, or the marginal log likelihood when marginal_likelihood is set.

Raises:

StateSpaceError: If \(F_t\) is not positive definite (code 2), the initialization is missing (3) or not stationary (6).

See also

loglike_concentrated: With the scale concentrated out.
filter: The full filter output.

loglike_concentrated()#

Evaluate the log likelihood with the scale concentrated out.

H, Q and \(P_*\) are taken relative to a scale \(\sigma^2\), which is replaced by its maximum likelihood estimate (DK §2.10.2, §7.3).

Returns:

llffloat: Concentrated log likelihood.
scalefloat: Estimate of \(\sigma^2\).

r2_diffuse(series=0)#

Compute the coefficient of determination against a random walk.

\(R^2_D = 1 - SSE / \sum_t (\Delta y_t - \overline{\Delta y})^2\) with SSE the sum of squared innovations after the diffuse periods (Harvey 1989; DK §7.4).

Parameters:

seriesint, optional: The series (0-based).

Returns:

float: \(R^2_D\); negative when the model forecasts worse than a random walk with drift.

simulate(rng=None, variates=None)#

Simulate observations and states from the model.

Parameters:

rngnumpy.random.Generator or int, optional: Source of the standard normal variates.
variatestuple of array_like, optional: The variates themselves, (u_init (m,), u_eps (p, n), u_eta (r, n)); overrides rng.

Returns:

yndarray, shape (p, n): Observations.
alphandarray, shape (m, n): States.
epsndarray, shape (p, n): Observation disturbances.
etandarray, shape (r, n): State disturbances.

Notes

\(\alpha_1 = a_1 + S(P_*) u_{init}\), so diffuse states start at their mean \(a_1\).

simulation_smoother(rng=None, variates=None)#

Draw states and disturbances from their distribution given the data.

Uses the mean-correction method of Durbin and Koopman (2002) (DK §4.9.2): simulate from the model, smooth the simulated data and shift by the smoothed values of the actual data.

Parameters:

rngnumpy.random.Generator or int, optional: Source of the standard normal variates.
variatestuple of array_like, optional: (u_init (m,), u_eps (p, n), u_eta (r, n)); overrides rng.

Returns:

alphandarray, shape (m, n): Draw of the states.
epsndarray, shape (p, n): Draw of the observation disturbances.
etandarray, shape (r, n): Draw of the state disturbances.

See also

djs_simulation_smoother: The de Jong-Shephard method (DK §4.9.3).

smooth(method='conventional')#

Run the filter and the state and disturbance smoother.

Parameters:

method{“conventional”, “sqrt”}, optional: "sqrt" uses the square root filter and smoother (DK §6.3).

Returns:

SmootherResults: Smoothed states and disturbances with their variances.

See also

fast_smoother, classical_smoother, two_filter_smoother: Alternative smoothers (DK §4.6).

Examples

>>> rep = ss.Representation(np.arange(5.0), k_states=1)
>>> rep["design"] = rep["transition"] = rep["selection"] = [[1.0]]
>>> rep["obs_cov"] = rep["state_cov"] = [[1.0]]
>>> rep.initialize_diffuse()
>>> res = rep.smooth()
>>> np.round(res.smoothed_state, 3)
array([[0.6, 1.2, 2. , 2.8, 3.4]])

smoothed_state_autocov()#

Compute the lag-one covariances of the smoothed state (DK §4.7).

Returns:

ndarray, shape (m, m, n - 1): Slice t is \(Cov(\alpha_{t+1}, \alpha_t | Y_n)\); NaN where the diffuse periods are involved.

smoothed_state_cov_between(t, j)#

Compute the covariance of two states given the data (DK §4.7).

Parameters:

t, jint: Periods (0-based), past the diffuse periods.

Returns:

ndarray, shape (m, m): \(Cov(\alpha_t, \alpha_j | Y_n)\).

smoothed_state_weights()#

Compute the weights of the data in the smoothed state (DK §4.8.3).

The smoothed state is linear in the data, the state intercepts and the initial mean: \(\hat\alpha_t = \sum_j W_{tj} (y_j - d_j) + \sum_j C_{tj} c_j + A_t a_1\).

Returns:

Wndarray, shape (m, p, n, n): W[:, :, t, j] is \(W_{tj}\).
Cndarray, shape (m, m, n, n): C[:, :, t, j] is \(C_{tj}\).
Andarray, shape (m, m, n): A[:, :, t] is \(A_t\).

Notes

Costs O(n) smoother runs.

steady_state()#

Compute the steady state of the Kalman filter (DK §2.11, §4.3.4).

The limit \(\bar P\) of the recursion for \(P_t\), by the structure-preserving doubling algorithm (Chu, Fan, Lin and Wang 2004), or by the recursion itself when H is singular.

Returns:

Pndarray, shape (m, m): \(\bar P\).
Fndarray, shape (p, p): \(\bar F = Z \bar P Z' + H\), the prediction error variance (DK §7.4).

Raises:

StateSpaceError: With code 4 if a system matrix varies over time, or 7 if the recursion does not converge.

Examples

The local level model with signal-to-noise ratio q has \(\bar P = x \sigma^2_\varepsilon\), \(x = (q + \sqrt{q^2 + 4q}) / 2\) (DK §2.11):

>>> rep = ss.Representation(np.zeros(10), k_states=1)
>>> rep["design"] = rep["transition"] = rep["selection"] = [[1.0]]
>>> rep["obs_cov"], rep["state_cov"] = [[1.0]], [[2.0]]
>>> P, F = rep.steady_state()
>>> q = 2.0
>>> float(P[0, 0]), float((q + np.sqrt(q**2 + 4 * q)) / 2)
(2.732050807568877, 2.732050807568877)

two_filter_smoother()#

Compute the smoothed state by the two-filter formula (DK §4.6.4).

A backward information filter combined with the forward filter. Not available with diffuse states.

Returns:

alphahatndarray, shape (m, n): Smoothed state.
Vndarray, shape (m, m, n): Its variance.

update_smoothed(alphahat, V, nobs_old)#

Update smoothed estimates for newly arrived observations (DK §4.4.5).

Parameters:

alphahatarray_like, shape (m, nobs_old): Smoothed state given the first nobs_old observations.
Varray_like, shape (m, m, nobs_old): Its variance.
nobs_oldint: Number of observations the estimates were based on.

Returns:

alphahatndarray, shape (m, n): Smoothed state given all n observations.
Vndarray, shape (m, m, n): Its variance.

whittle_smoother()#

Compute the smoothed state from the Whittle relation (DK §4.6.3).

Returns:

ndarray, shape (m, n): Smoothed state.

Notes

The recursion runs backwards through \(T_t^{-1}\) and loses accuracy with the length of the series, as DK note. It serves to illustrate the relation; use smooth in practice.