geossm.stmodel.LRStateSpaceResults

class geossm.stmodel.LRStateSpaceResults(model=None, params=None, nstats=None, options=None, y_hat_list=None, Sigma_y_hat_list=None, block_p=None, points_hat=None, timestamps_hat=None, crs_hat=None, **kwargs)

Bases: StateSpaceResults

Results of fitting a geossm.stmodel.LRStateSpaceModel.

Subclass of geossm.ssm.StateSpaceResults returned by LRStateSpaceModel.fit(), carrying the estimated ModelParams (.params), the EM iteration history (.nstats/.llf_path), and the in-sample fitted values/smoothed states, in addition to everything the base class already provides (residuals, mse/rmse, conf_int_y, …). It also adds:

• Parameter inference: compute_cov_params (observed-information Hessian, via jax.hessian on the model’s log-likelihood), bse, tvalues, pvalues, conf_int, aic/bic.

• Out-of-sample prediction: predict (thin wrapper around LRStateSpaceModel.predict, storing y_pred_list/ Sigma_y_pred_list on this object).

• Back-transformation of predictions/fitted values to the response’s original scale via the delta method (back_transform), when the model formula applies a transform to the response (e.g. np.log).

• Export to geopandas.GeoDataFrame (to_geo).

• A statsmodels-style textual report (summary).

Typical usage, continuing from LRStateSpaceModel.fit():

results = model.fit()
results.compute_cov_params()      # optional, needed for bse/t/p/CI
print(results.summary())
results = results.predict(new_df) # out-of-sample prediction
results = results.back_transform(g_inv=jnp.exp)  # if y = np.log(...)
geo = results.to_geo()

Methods

as_dict	Return a plain dict snapshot of the main results.
back_transform	Map y_hat_list/y_pred_list (and their model-implied covariance) back to the response's original scale, via the second-order delta method (see _delta_method), applied per response variable, and store the results as y_hat_back_list/Sigma_y_hat_back_list and -- if .predict() has been run -- y_pred_back_list/ Sigma_y_pred_back_list.
compute_aic	Akaike Information Criterion: 2 * k - 2 * llf, with k = n_params (free parameters) and llf the fitted log-likelihood.
compute_bic	Bayesian Information Criterion: log(nobs) * k - 2 * llf, with k = n_params (free parameters), nobs the number of observed values, and llf the fitted log-likelihood.
compute_cov_params	Compute the asymptotic covariance matrix of the free parameters from the observed-information Hessian of the log-likelihood (-inverse(Hessian), symmetrised), and cache it.
conf_int	Confidence intervals for the free parameters.
conf_int_state	Compute pointwise Gaussian confidence intervals for the latent state.
conf_int_y	Compute pointwise Gaussian confidence/prediction intervals for the predicted observations y_hat.
coverage_probability	Empirical coverage probability of the conf_int_y prediction intervals.
diagnostics	Compute standard residual-normality/autocorrelation diagnostics.
generate_summary	Build the left/right key-value rows used by summary()'s header table: the base class's rows (model info, fit quality) plus the EM iteration count/runtime and AIC/BIC, and -- if .predict() has been run -- a small summary of the out-of-sample predictions (see _pred_summary_stats).
mse	Mean squared error of the residuals.
predict	Compute out-of-sample predictions based on smoothed states and model parameters.
rmse	Root mean squared error of the residuals.
summary	Build a statsmodels-style textual report of the fitted model: the header table from generate_summary (model info, EM/fit diagnostics, AIC/BIC, and prediction summary if available), followed by one parameter table per group -- fixed effects (beta), measurement-equation parameters (s2e, A), and state-equation parameters (ks, f) -- each with point estimate, standard error, t-value, p-value, and alpha-level confidence interval.
to_geo	Package the in-sample fitted values and, if .predict() has been run, the out-of-sample prediction into GeoDataFrames, one row per (point, timestamp), ready to be exported (e.g. .to_file("out.shp")).
to_numpy	Convert array-like attributes to numpy arrays in the results' dtype (the model's dtype, by default), so results stay numerically consistent with the precision the model was fit in.
update	Return a new StateSpaceResults with some fields replaced.

Attributes

aic	Akaike Information Criterion, compute_aic().
backend	JAX compute device any further computation on these results runs on (defaults to the originating model's geossm.ssm.statespace.StateSpaceModel.backend).
bic	Bayesian Information Criterion, compute_bic().
bse	Standard errors of the free parameters, reshaped back into a ModelParams structure (one array per parameter, stored in each Param.bse) via _vector_to_bse_params, mirroring the shape of self.params.
bse_vector	Standard errors of the free parameters, as a flat 1D array (sqrt of the diagonal of compute_cov_params(), clipped at 0 to guard against small negative values from numerical noise), in the same flat order as theta_hat.
df_resid	Residual degrees of freedom, nobs - n_params (number of observed values minus the number of free parameters).
dtype	Floating-point precision results are stored in (defaults to the originating model's dtype).
n_params	Total number of free (not fixed) scalar parameters across beta, s2e, f, A, ks, x0, Sigma0, used as the degrees of freedom / penalty term k in df_resid, compute_aic, and compute_bic.
pvalues	Two-sided p-values for the free parameters (Student-t with df_resid degrees of freedom), in the same flat order as theta_hat.
residuals	Observation residuals, or None if unavailable.
residuals_list	Per-variable view of the base class's stacked residuals (y_obs - y_hat), split along block_p -- the in-sample counterpart of y_hat_list/Sigma_y_hat_list.
theta_hat	Point estimates of the free (not fixed) parameters, flattened into a single 1D array in the order produced by _pack_params (x0/Sigma0 always treated as fixed here, see _inference_params) -- the same layout as bse_vector, tvalues, pvalues, and conf_int.
tvalues	t-statistics (theta_hat / bse_vector) for the free parameters, in the same flat order as theta_hat.

property n_params

Total number of free (not fixed) scalar parameters across beta, s2e, f, A, ks, x0, Sigma0, used as the degrees of freedom / penalty term k in df_resid, compute_aic, and compute_bic.

Type:

int

compute_cov_params()

Compute the asymptotic covariance matrix of the free parameters from the observed-information Hessian of the log-likelihood (-inverse(Hessian), symmetrised), and cache it.

This drives bse/bse_vector, tvalues, pvalues, and conf_int, and its presence (together with summary()) decides whether summary() shows real standard errors or NaN placeholders (see _nan_bse_params). Computing the Hessian (_compute_hessian, via jax.hessian) can be relatively slow for models with many free parameters, so it is not run automatically by fit(); call this explicitly when standard errors are needed.

Returns:

jax.numpy.ndarray – The covariance matrix of the free (not fixed) parameters, in the flat order given by _pack_params/theta_hat. Cached on self._cov_params after the first call.

property df_resid

Residual degrees of freedom, nobs - n_params (number of observed values minus the number of free parameters). Used as the Student-t degrees of freedom in _stats_from_arrays (p-values and confidence intervals).

Type:

int

property bse_vector

Standard errors of the free parameters, as a flat 1D array (sqrt of the diagonal of compute_cov_params(), clipped at 0 to guard against small negative values from numerical noise), in the same flat order as theta_hat. Triggers the Hessian computation on first access if not already cached.

Type:

ndarray

property bse

Standard errors of the free parameters, reshaped back into a ModelParams structure (one array per parameter, stored in each Param.bse) via _vector_to_bse_params, mirroring the shape of self.params. Fixed parameters get bse=None.

Type:

ModelParams

property theta_hat

Point estimates of the free (not fixed) parameters, flattened into a single 1D array in the order produced by _pack_params (x0/Sigma0 always treated as fixed here, see _inference_params) – the same layout as bse_vector, tvalues, pvalues, and conf_int.

Type:

ndarray

property tvalues

t-statistics (theta_hat / bse_vector) for the free parameters, in the same flat order as theta_hat.

Type:

ndarray

property pvalues

Two-sided p-values for the free parameters (Student-t with df_resid degrees of freedom), in the same flat order as theta_hat.

Type:

ndarray

conf_int(alpha=0.05)

Confidence intervals for the free parameters.

Parameters:

alpha (float, default 0.05) – Significance level; the returned interval has coverage 1 - alpha (e.g. alpha=0.05 gives a 95% confidence interval).

Returns:

ndarray, shape (n_free_params, 2) – Lower/upper confidence bounds, in the same flat parameter order as theta_hat.

property aic

Akaike Information Criterion, compute_aic().

Type:

float

property bic

Bayesian Information Criterion, compute_bic().

Type:

float

compute_aic()

Akaike Information Criterion: 2 * k - 2 * llf, with k = n_params (free parameters) and llf the fitted log-likelihood.

Raises:

AttributeError – If llf is not set (only available when nstats was provided to __init__, i.e. after LRStateSpaceModel.fit()).

compute_bic()

Bayesian Information Criterion: log(nobs) * k - 2 * llf, with k = n_params (free parameters), nobs the number of observed values, and llf the fitted log-likelihood.

Raises:

AttributeError – If llf is not set (only available when nstats was provided to __init__, i.e. after LRStateSpaceModel.fit()).

predict(df, verbose=True)

Compute out-of-sample predictions based on smoothed states and model parameters. self.model.predict stores them directly on this results object (points_pred, y_pred_list, Sigma_y_pred_list, tdelta_pred, timestamps_pred, crs_pred – one entry per response variable) and returns self, so predictions travel with the fitted model and can be reused by other methods (e.g. .to_geo(), plotting, summaries) without re-running prediction.

Parameters:

• df (geopandas.GeoDataFrame) – New locations/times to predict at; see LRStateSpaceModel.predict for the exact requirements.

• verbose (bool, default True) – Verbosity for the underlying grid/design-matrix construction.

Returns:

LRStateSpaceResults – self, with the prediction attributes populated.

property residuals_list

Per-variable view of the base class’s stacked residuals (y_obs - y_hat), split along block_p – the in-sample counterpart of y_hat_list/Sigma_y_hat_list.

to_geo()

Package the in-sample fitted values and, if .predict() has been run, the out-of-sample prediction into GeoDataFrames, one row per (point, timestamp), ready to be exported (e.g. .to_file(“out.shp”)).

Returns a dict with two keys:

• “hat”: GeoDataFrame over the training grid, from y_hat_list/Sigma_y_hat_list (plus y_hat_back_<var>/ std_hat_back_<var> if .back_transform() has been run). Always present.

• “pred”: GeoDataFrame over the prediction grid, from y_pred_list/Sigma_y_pred_list (plus y_pred_back_<var>/ std_pred_back_<var> if .back_transform() has been run), or None if .predict() hasn’t been run yet.

These are two separate GeoDataFrames, not one, because the training and prediction grids generally have a different number of points/timestamps.

Columns: point_id, timestamp, then y_<hat|pred>_<var>/ std_<hat|pred>_<var> for every response variable.

results = results.predict(grid, verbose=True) geo = results.to_geo() geo[“hat”].to_file(“fitted.shp”) geo[“pred”].to_file(“predictions.shp”)

as_dict()

Return a plain dict snapshot of the main results.

Ensures all array attributes are NumPy arrays first (via to_numpy). Useful for serialization or for passing results to code that should not depend on the StateSpaceResults class.

Returns:

dict – Keys: y_obs, y_hat, Sigma_y_hat, x_filtered, P_filtered, x_smoothed, P_smoothed, llf, mse (global), rmse (global), diagnostics (see diagnostics), S11, S10, S00.

Return type:

Dict[str, Any]

back_transform(g_inv)

Map y_hat_list/y_pred_list (and their model-implied covariance) back to the response’s original scale, via the second-order delta method (see _delta_method), applied per response variable, and store the results as y_hat_back_list/Sigma_y_hat_back_list and – if .predict() has been run – y_pred_back_list/ Sigma_y_pred_back_list.

This relies on the transformed response being asymptotically Normal (the model’s own assumption), so the delta method’s local, second-order Taylor expansion around the mean is a reasonable approximation – but it is not the only option. Alternatives, roughly in order of how much they trade simplicity for accuracy:

• Naive plug-in g_inv(y_hat): what the first-order term alone gives; biased whenever g_inv is curved (Jensen’s inequality), which is exactly what the second-order correction here fixes.

• Closed-form formulas for a specific g_inv, when known – e.g. the lognormal mean exp(mu + sigma^2/2) for log, which is the delta method’s answer and the exact one for that case.

• Duan’s smearing estimator: replace h(mu) by the empirical average of h(mu + residual_i) over the in-sample residuals. Distribution-free (no normality needed) but needs those residuals kept around, and is usually applied to the mean only.

• Gauss-Hermite quadrature / Monte Carlo against Normal(mu, Var): numerically integrates h under the same normality assumption used here, so it stays exact as curvature or Var grows large (where a second-order Taylor expansion starts to break down), at the cost of a handful of extra evaluations of h per point.

The delta method is the right default when Var is small relative to h’s curvature (the usual case for a well-identified model); if predictions look off for locations/times with large predictive variance, quadrature is the natural drop-in replacement since it reuses the same mu/Sigma this method already computes.

Parameters:

g_inv (Callable[[float], float]) – The inverse of the response transform in the model formula (e.g. jnp.exp for a np.log(y) response), as a scalar -> scalar JAX-traceable function.

property backend

JAX compute device any further computation on these results runs on (defaults to the originating model’s geossm.ssm.statespace.StateSpaceModel.backend).

Type:

jax.Device

conf_int_state(which='smoothed', alpha=0.05)

Compute pointwise Gaussian confidence intervals for the latent state.

For each time step, uses the marginal state mean/variance (\(x_t\), \(\mathrm{diag}(P_t)\)) and a two-sided Normal quantile: \(x_t \pm z_{1-\alpha/2} \sqrt{\mathrm{diag}(P_t)}\). Does not account for cross-state covariances (only the marginal per-component variance is used).

Parameters:

• which ({'smoothed', 'filtered'}, default 'smoothed') – Whether to use the smoothed (x_smoothed/P_smoothed) or filtered (x_filtered/P_filtered) state estimates.

• alpha (float, default 0.05) – Significance level; the interval has nominal coverage 1 - alpha.

Returns:

lower, upper (numpy.ndarray, shape (q, T+1)) – Lower and upper confidence bounds for the state mean.

Raises:

ValueError – If which is invalid, or the required state means/ covariances are not available on this results object.

Return type:

Tuple[ndarray, ndarray]

conf_int_y(alpha=0.05, prediction=False)

Compute pointwise Gaussian confidence/prediction intervals for the predicted observations y_hat.

For each time step, uses \(\hat y_t \pm z_{1-\alpha/2} \sqrt{\mathrm{diag}(\Sigma_{\hat y, t})}\), where \(\Sigma_{\hat y, t}\) is Sigma_y_hat (the state-driven predictive covariance, excluding observation noise) plus, when prediction=True, the observation noise covariance R read from self.model - turning a confidence interval for the mean response into a prediction interval for a new observation.

Parameters:

• alpha (float, default 0.05) – Significance level; the interval has nominal coverage 1 - alpha.

• prediction (bool, default False) – If True, add the model’s observation noise covariance R to Sigma_y_hat before computing the interval (prediction interval); if False, use Sigma_y_hat alone (confidence interval for the mean).

Returns:

lower, upper (numpy.ndarray, shape (p, T)) – Lower and upper interval bounds for y_hat.

Raises:

ValueError – If y_hat or Sigma_y_hat is not available on this results object.

Return type:

Tuple[ndarray, ndarray]

coverage_probability(alpha=0.05, which='global')

Empirical coverage probability of the conf_int_y prediction intervals.

Computes the fraction of observed values y_obs that fall within the 1 - alpha prediction interval from conf_int_y(alpha, prediction=True); for a well-calibrated model this should be close to 1 - alpha.

Parameters:

• alpha (float, default 0.05) – Significance level used for the prediction interval.

• which ({'global', 'space', 'time'}, default 'global') – Aggregation axis: 'global' averages over all entries; 'space' averages over time (shape (p,)); 'time' averages over space (shape (T,)).

Returns:

float or numpy.ndarray – Empirical coverage probability in [0, 1], aggregated as requested.

diagnostics()

Compute standard residual-normality/autocorrelation diagnostics.

Flattens the (non-missing) residuals and computes:

• Jarque-Bera test statistic and p-value (jb/jb_pvalue) for normality, based on sample skewness and kurtosis.

• Omnibus normality test statistic and p-value (omni/omni_pvalue, statsmodels’ omni_normtest).

• Sample skewness and excess kurtosis (skew/kurtosis).

• Durbin-Watson statistic (dw) for first-order autocorrelation (values near 2 indicate little autocorrelation).

Returns:

dict of str to float – Keys jb, jb_pvalue, omni, omni_pvalue, dw, skew, kurtosis.

Return type:

Dict[str, float]

property dtype

Floating-point precision results are stored in (defaults to the originating model’s dtype).

Type:

numpy.dtype

mse(which='global')

Mean squared error of the residuals.

Parameters:

which ({'global', 'space', 'time'}, default 'global') – Aggregation axis: 'global' averages over all (non-missing) entries; 'space' averages over time, returning one value per observed location/variable, shape (p,); 'time' averages over space, returning one value per time step, shape (T,). NaN residuals (missing observations) are ignored.

Returns:

float or numpy.ndarray – Mean squared error, aggregated as requested.

Raises:

ValueError – If which is not one of 'global', 'space', 'time'.

Return type:

float

property residuals: ndarray | None

Observation residuals, or None if unavailable.

Returns:

numpy.ndarray or None – y_obs - y_hat, shape (p, T) (NaN where y_obs is missing), or None if either y_obs or y_hat is unavailable.

rmse(which='global')

Root mean squared error of the residuals.

Parameters:

which ({'global', 'space', 'time'}, default 'global') – Aggregation axis; forwarded to mse.

Returns:

float or numpy.ndarray – Square root of mse, aggregated as requested.

Return type:

float

to_numpy()

Convert array-like attributes to numpy arrays in the results’ dtype (the model’s dtype, by default), so results stay numerically consistent with the precision the model was fit in.

update(**kwargs)

Return a new StateSpaceResults with some fields replaced.

Copies every current public (non-underscore-prefixed) attribute except today, overrides the ones given in kwargs, and constructs a new instance from the result. Does not mutate self.

Parameters:

**kwargs – Attribute name/value pairs to override; each key must already be a public attribute of this instance.

Returns:

StateSpaceResults – A new instance with the requested fields updated.

Raises:

AttributeError – If a key in kwargs is not an existing public attribute.

Return type:

StateSpaceResults

generate_summary()

Build the left/right key-value rows used by summary()’s header table: the base class’s rows (model info, fit quality) plus the EM iteration count/runtime and AIC/BIC, and – if .predict() has been run – a small summary of the out-of-sample predictions (see _pred_summary_stats).

Returns:

gen_top_left, gen_top_right (list of (str, list)) – Two lists of (label, [value]) rows, of equal length.

summary(alpha=0.05)

Build a statsmodels-style textual report of the fitted model: the header table from generate_summary (model info, EM/fit diagnostics, AIC/BIC, and prediction summary if available), followed by one parameter table per group – fixed effects (beta), measurement-equation parameters (s2e, A), and state-equation parameters (ks, f) – each with point estimate, standard error, t-value, p-value, and alpha-level confidence interval.

If compute_cov_params() has not been called yet, standard errors/t-values/p-values/confidence intervals are shown as NaN placeholders (see _nan_bse_params) rather than triggering the (potentially slow) Hessian computation implicitly.

Parameters:

alpha (float, default 0.05) – Significance level for the parameter tables’ confidence intervals (see conf_int).

Returns:

statsmodels.iolib.summary.Summary – Printable summary object (str(…)/print(…)).