BayesianVARMAX#

class pymc_experimental.statespace.models.BayesianVARMAX(order: Tuple[int, int], endog_names: list[str] | None = None, k_endog: int | None = None, stationary_initialization: bool = False, filter_type: str = 'standard', measurement_error: bool = False, verbose=True)[source]#

Vector AutoRegressive Moving Average with eXogenous Regressors

Parameters:

order (tuple of (int, int)) – Number of autoregressive (AR) and moving average (MA) terms to include in the model. All terms up to the specified order are included. For restricted models, set zeros directly on the priors.
endog_names (list of str, optional) –
Names of the endogenous variables being modeled. Used to generate names for the state and shock coords. If None, the state names will simply be numbered.

Exactly one of either endog_names or k_endog must be specified.
k_endog (int, optional) –
Number of endogenous states to be modeled.

Exactly one of either endog_names or k_endog must be specified.
stationary_initialization (bool, default False) –
If true, the initial state and initial state covariance will not be assigned priors. Instead, their steady state values will be used. If False, the user is responsible for setting priors on the initial state and initial covariance.

..warning :: This option is very sensitive to the priors placed on the AR and MA parameters. If the model dynamics
for a given sample are not stationary, sampling will fail with a “covariance is not positive semi-definite” error.
filter_type (str, default "standard") – The type of Kalman Filter to use. Options are “standard”, “single”, “univariate”, “steady_state”, and “cholesky”. See the docs for kalman filters for more details.
state_structure (str, default "fast") –
How to represent the state-space system. When “interpretable”, each element of the state vector will have a precise meaning as either lagged data, innovations, or lagged innovations. This comes at the cost of a larger state vector, which may hurt performance.

When “fast”, states are combined to minimize the dimension of the state vector, but lags and innovations are mixed together as a result. Only the first state (the modeled timeseries) will have an obvious interpretation in this case.
measurement_error (bool, default True) – If true, a measurement error term is added to the model.
verbose (bool, default True) – If true, a message will be logged to the terminal explaining the variable names, dimensions, and supports.

Notes

The VARMA model is a multivariate extension of the SARIMAX model. Given a set of timeseries \(\{x_t\}_{t=0}^T\), with \(x_t = \begin{bmatrix} x_{1,t} & x_{2,t} & \cdots & x_{k,t} \end{bmatrix}^T\), a VARMA models each series as a function of the histories of all series. Specifically, denoting the AR-MA order as (p, q), a VARMA can be written:

\[x_t = A_1 x_{t-1} + A_2 x_{t-2} + \cdots + A_p x_{t-p} + B_1 \varepsilon_{t-1} + \cdots + B_q \varepsilon_{t-q} + \varepsilon_t\]

Where \(\varepsilon_t = \begin{bmatrix} \varepsilon_{1,t} & \varepsilon_{2,t} & \cdots & \varepsilon_{k,t}\end{bmatrix}^T \sim N(0, \Sigma)\) is a vector of i.i.d stochastic innovations or shocks that drive intertemporal variation in the data. Matrices \(A_i, B_i\) are \(k \times k\) coefficient matrices:

\[\begin{split}A_i = \begin{bmatrix} \rho_{1,i,1} & \rho_{1,i,2} & \cdots & \rho_{1,i,k} \\ \rho_{2,i,1} & \rho_{2,i,2} & \cdots & \rho_{2,i,k} \\ \vdots & \vdots & \cdots & \vdots \\ \rho{k,i,1} & \rho_{k,i,2} & \cdots & rho_{k,i,k} \end{bmatrix}\end{split}\]

Internally, this representation is not used. Instead, the vectors \(x_t, x_{t-1}, \cdots, x_{t-p}, \varepsilon_{t-1}, \cdots, \varepsilon_{t-q}\) are concatenated into a single column vector of length k * (p+q), while the coefficients matrices are likewise concatenated into a single coefficient matrix, \(T\).

As the dimensionality of the VARMA system increases – either because there are a large number of timeseries included in the analysis, or because the order is large – the probability of sampling a stationary matrix \(T\) goes to zero. This has two implications for applied work. First, a non-stationary system will exhibit explosive behavior, potentially rending impulse response functions and long-term forecasts useless. Secondly, it is not possible to do stationary initialization. Stationary initialization significantly speeds up sampling, and should be preferred when possible.

Examples

The following code snippet estimates a VARMA(1, 1):

import pymc_experimental.statespace as pmss
import pymc as pm

# Create VAR Statespace Model
bvar_mod = pmss.BayesianVARMAX(endog_names=data.columns, order=(2, 0),
                               stationary_initialization=False, measurement_error=False,
                               filter_type="standard", verbose=True)

# Unpack dims and coords
x0_dims, P0_dims, state_cov_dims, ar_dims = bvar_mod.param_dims.values()
coords = bvar_mod.coords

# Estimate PyMC model
with pm.Model(coords=coords) as var_mod:
    x0 = pm.Normal("x0", dims=x0_dims)
    P0_diag = pm.Gamma("P0_diag", alpha=2, beta=1, size=data.shape[1] * 2, dims=P0_dims[0])
    P0 = pm.Deterministic("P0", pt.diag(P0_diag), dims=P0_dims)

    state_chol, _, _ = pm.LKJCholeskyCov(
        "state_chol", eta=1, n=bvar_mod.k_posdef, sd_dist=pm.Exponential.dist(lam=1)
    )

    ar_params = pm.Normal("ar_params", mu=0, sigma=1, dims=ar_dims)
    state_cov = pm.Deterministic("state_cov", state_chol @ state_chol.T, dims=state_cov_dims)

    bvar_mod.build_statespace_graph(data, mode="JAX")
    idata = pm.sample(nuts_sampler="numpyro")

__init__(order: Tuple[int, int], endog_names: list[str] | None = None, k_endog: int | None = None, stationary_initialization: bool = False, filter_type: str = 'standard', measurement_error: bool = False, verbose=True)[source]#

Methods

`__init__`(order[, endog_names, k_endog, ...])
`add_default_priors`()	Add default priors to the active PyMC model context
`add_exogenous`(exog)	Add an exogenous process to the statespace model
`build_statespace_graph`(data[, ...])	Given a parameter vector theta, constructs the full computational graph describing the state space model and the associated log probability of the data.
`forecast`(idata, start[, periods, end, ...])	Generate forecasts of state space model trajectories into the future.
`impulse_response_function`(idata[, n_steps, ...])	Generate impulse response functions (IRF) from state space model dynamics.
`make_and_register_data`(name, shape[, dtype])	Helper function to create a pytensor symbolic variable and register it in the _name_to_data dictionary
`make_and_register_variable`(name, shape[, dtype])	Helper function to create a pytensor symbolic variable and register it in the _name_to_variable dictionary
`make_symbolic_graph`()	The purpose of the make_symbolic_graph function is to hide tedious parameter allocations from the user.
`sample_conditional_posterior`(idata[, ...])	Sample from the conditional posterior; that is, given parameter draws from the posterior distribution, compute Kalman filtered trajectories.
`sample_conditional_prior`(idata[, random_seed])	Sample from the conditional prior; that is, given parameter draws from the prior distribution, compute Kalman filtered trajectories.
`sample_unconditional_posterior`(idata[, ...])	Draw unconditional sample trajectories according to state space dynamics, using random samples from the posterior distribution over model parameters.
`sample_unconditional_prior`(idata[, steps, ...])	Draw unconditional sample trajectories according to state space dynamics, using random samples from the prior distribution over model parameters.
`unpack_statespace`()	Helper function to quickly obtain all statespace matrices in the standard order.

Attributes

`coords`	PyMC model coordinates
`data_info`	Information about Data variables that need to be declared in the PyMC model block.
`data_names`	Names of data variables expected by the model.
`default_priors`	Dictionary of parameter names and callable functions to construct default priors for the model
`observed_states`	A k_endog length list of strings, associated with the model's observed states
`param_dims`	Dictionary of named dimensions for each model parameter
`param_info`	Information about parameters needed to declare priors
`param_names`	Names of model parameters
`shock_names`	A k_posdef length list of strings, associated with the model's shock processes
`state_names`	A k_states length list of strings, associated with the model's hidden states