Observation-driven recursions

.. py:module:: pymc_forecast.ssoe

Observation-driven recursions with a single source of error.

The deterministic training filter and the generative forecast share the same mean and update functions. Only future errors are random variables; they are absent during fitting and drawn fresh during posterior replay. Inspired by numpyro_forecast.models.ssoe (Apache-2.0), using PyTensor and named coords.

.. py:class:: SSOEResult(mu, mu_future, y_future, dims) :module: pymc_forecast.ssoe

Bases: :py:class:object

Outputs of :func:ssoe, before registering an observation likelihood.

mu contains in-sample one-step-ahead means. mu_future contains forecast means conditional on the preceding simulated observations; y_future contains those means plus their current observation errors. Both future tensors have length zero during training. Time is first in these symbolic tensors; dims records the remaining named dimensions. Use mu_future to exclude the current observation error, remembering that earlier future errors still affect the state.

Each parameter is available as an attribute of the same name.

:param mu: In-sample one-step-ahead means, shape (t_obs, *row). :type mu: :py:class:pytensor.tensor.TensorVariable :param mu_future: Future means, shape (future, *row) (length zero during training). :type mu_future: :py:class:pytensor.tensor.TensorVariable :param y_future: Future simulated observations mu_future plus the current errors (length zero during training). :type y_future: :py:class:pytensor.tensor.TensorVariable :param dims: Non-time dims of the rows, in the order used by the tensors. :type dims: :py:class:tuple of :py:class:str

.. py:attribute:: SSOEResult.mu :module: pymc_forecast.ssoe

.. py:attribute:: SSOEResult.mu_future :module: pymc_forecast.ssoe

.. py:attribute:: SSOEResult.y_future :module: pymc_forecast.ssoe

.. py:attribute:: SSOEResult.dims :module: pymc_forecast.ssoe

.. py:function:: ssoe(h, name, y, init, mean, update, noise, xs=None, *, params=(), dims=None) :module: pymc_forecast.ssoe

Filter observed values, then simulate a recursive forecast.

mean(state, x, *params) and update(state, y, error, x, *params) take params because PyTensor scan cannot close over random variables. noise is not scanned: it is a .dist() or a Prior registered only as {name}_future. Pass random coefficients in params. This helper does not auto-detect closed-over RVs.

:param h: The horizon of the current model build. :type h: :py:class:pymc_forecast.model.Horizon :param name: Base name of the future error variable. Only f"{name}_future" is registered, and only when forecasting. :type name: :py:class:str :param y: Labeled driving history. None uses h.data. Must cover exactly the training window (time coords equal to h.time), with finite values and non-time coords matching the model coords. A transformed history can be supplied to compose multiple recursion channels. Prior-only builds need an explicit driving history; this helper does not generate an in-sample history. :type y: :py:class:xarray.DataArray or :py:obj:None :param init: Initial state: one tensor-like or a nonempty tuple of them (only a tuple means multiple states; each is cast to floatX). States may have different shapes (e.g. scalar level and vector seasonality). :type init: :py:class:pytensor.tensor.TensorVariable or :py:class:tuple of :py:class:pytensor.tensor.TensorVariable :param mean: (state, x_t, *params) -> mu_t. Returns the one-step-ahead mean, shaped like one row of y. x_t is None without xs. :type mean: :py:class:callable :param update: (state, y_t, eps_t, x_t, *params) -> state. Returns the next state with the same structure and shapes as init. In-sample, eps_t = y_t - mu_t; in the future, eps_t is freshly drawn and y_t = mu_t + eps_t. Both callbacks must be deterministic. :type update: :py:class:callable :param noise: Unnamed .dist() or pymc-extras Prior for independent, zero-centered per-step future errors, expanded and registered as f"{name}_future" only. An RVFactory or other callable is rejected. Parameters broadcast against ("time_future", *dims), so a per-series scale works. Errors may be correlated across the observation dimensions, using e.g. pm.MvNormal.dist. In-sample errors are residuals, not random variables.

             A learned scale must be a model variable shared with the observation:
             create ``sigma = pm.HalfNormal("sigma", 1)`` in the model body, use it
             in the observation, and pass ``pm.Normal.dist(0, sigma)``. Hence a
             ``Prior`` with a hyper-prior parameter (which would exist only on
             forecasting builds and be drawn from its prior), a model variable
             instead of a ``.dist()``, and a ``.dist()`` depending on unnamed random
             variables are rejected, on training builds too. A Normal/StudentT
             ``Prior`` must leave ``mu`` unset or zero, and a ``pm.Normal.dist`` /
             ``pm.StudentT.dist`` must have a constant zero location
             (``pm.Normal.dist(sigma)`` binds ``sigma`` as ``mu``). The locations of
             other dists are not checked.

:type noise: :py:class:pymc_extras.prior.Prior or :py:class:pytensor.tensor.TensorVariable :param xs: Labeled inputs spanning the full horizon (None means no inputs). The time dimension is selected by name and its coordinates are checked; non-time dims are checked only when they are model coords. Values must be finite. Future inputs must be known covariates or explicit scenarios: never derive future update gates from held-out observations. Extra rows are ignored during a shorter training build. :type xs: :py:class:xarray.DataArray, optional :param params: Tensor parameters passed explicitly to both callbacks. Required for random coefficients: scan cannot close over RVs, and closed-over RVs are not detected. :type params: :py:class:sequence, default () :param dims: Non-time observation dimensions; None infers them from y. Labeled data are transposed into this order before entering the scan. :type dims: :py:class:tuple of :py:class:str, optional

:returns: :py:class:SSOEResult – In-sample means, future means and future samples. The caller registers obs against result.mu and forecast as a Deterministic of result.y_future. Register mu and mu_future Deterministics to include the means in the standard prediction outputs. Do not add another observation draw to y_future: it already includes noise.

:raises pymc_forecast.exceptions.HorizonError: If there is no observed history (y and h.data are both None, or h.t_obs == 0); if noise is invalid (see above; checked on every build); or, when forecasting, if a dim in ("time_future", *dims) is not a model coord or a .dist() noise has more axes than dims. :raises pymc_forecast.exceptions.AlignmentError: If y or xs is not a DataArray with a "time" dim; if the time coords of y differ from h.time, or xs is shorter than the horizon or its time coords differ; if dims does not name each non-time dim of y exactly once or contains "time_future"; or if a dim’s coords do not match the model coords. :raises ValueError: If y or xs has non-finite values; if init is an empty tuple; or if mean/update violate their contract (mean with the wrong number of dims, update changing the state structure, count or number of dims, or the callbacks creating random updates or model variables).

.. rubric:: Notes

This is an observation-driven filter, not a latent Markov process. For sampled hidden states use :func:~pymc_forecast.markov.markov_series; for linear-Gaussian hidden states consider the statespace backend (:class:~pymc_forecast.statespace.StatespaceForecaster).

A minimal ARMA model

The state contains the previous observation and its one-step prediction error. The same callbacks filter training observations and simulate future observations.

import pymc as pm
from pymc_forecast import Horizon, ssoe


def arma(covariates, data=None):
    h = Horizon.from_data(covariates, data)
    phi = pm.Uniform("phi", -0.95, 0.95)
    theta = pm.Uniform("theta", -0.95, 0.95)
    sigma = pm.HalfNormal("sigma", 1)

    result = ssoe(
        h,
        "eps",
        None,
        (0.0, 0.0),
        lambda state, x, phi, theta: phi * state[0] + theta * state[1],
        lambda state, y, error, x, phi, theta: (y, error),
        pm.Normal.dist(0, sigma),
        params=(phi, theta),
    )
    pm.Normal("obs", result.mu, sigma, observed=h.data.values, dims="time")
    pm.Deterministic("mu", result.mu, dims="time")
    if h.future:
        pm.Deterministic("mu_future", result.mu_future, dims="time_future")
        pm.Deterministic("forecast", result.y_future, dims="time_future")

The helper registers only eps_future, and only during forecasting. The caller owns the observation distribution and can compose several recursion channels. Do not pass y_future to predict with another noisy observation distribution: y_future already contains the observation error used to update the future state. The explicit Deterministics above use the usual prediction schema and are returned by forecast and predict_in_sample.

mu is a one-step fitted mean conditioned on the actual observed history. In-sample posterior prediction samples observations around these fitted means; it does not simulate a new training history. mu_future conditions on previously simulated future observations and excludes only the current step’s error. It is not an unconditional multi-step mean with all future errors integrated out.

Labeled inputs and update gates

y defaults to h.data; a transformed DataArray can be passed explicitly. Its time coordinates must match the observed window exactly. Non-time dimensions are inferred by name, checked against model coordinates, and preserved in the result. xs is a DataArray covering the full horizon; its time dimension can appear anywhere. Each callback receives one time slice in its other dimension order. Pass model parameters via params and keep both callbacks deterministic.

For unavailable or missing observations, supply a finite placeholder and use an explicit gate in xs to prevent state updates. The likelihood must also exclude those observations; an update gate alone does not mask a likelihood. Future gates must come from known inputs or an explicit scenario, never held-out observations. The helper does not infer a missing-data or censoring mechanism.

For sampled hidden states, use markov_series. For linear-Gaussian hidden states, the pymc-extras statespace backend can marginalize the state path. See the ARMA, Holt-Winters, and censored demand examples for observation-driven recursions. The censored example is the update-gate pattern above: availability keeps a stockout zero out of the lag, and the likelihood contributes survival mass at a known cap instead of a density.