class_docstr: Vector Autoregressive Moving Average with eXogenous regressors model
  
  Parameters
  ----------
  endog : array_like
      The observed time-series process :math:`y`, , shaped nobs x k_endog.
  exog : array_like, optional
      Array of exogenous regressors, shaped nobs x k.
  order : iterable
      The (p,q) order of the model for the number of AR and MA parameters to
      use.
  trend : str{'n','c','t','ct'} or iterable, optional
      Parameter controlling the deterministic trend polynomial :math:`A(t)`.
      Can be specified as a string where 'c' indicates a constant (i.e. a
      degree zero component of the trend polynomial), 't' indicates a
      linear trend with time, and 'ct' is both. Can also be specified as an
      iterable defining the non-zero polynomial exponents to include, in
      increasing order. For example, `[1,1,0,1]` denotes
      :math:`a + bt + ct^3`. Default is a constant trend component.
  error_cov_type : {'diagonal', 'unstructured'}, optional
      The structure of the covariance matrix of the error term, where
      "unstructured" puts no restrictions on the matrix and "diagonal"
      requires it to be a diagonal matrix (uncorrelated errors). Default is
      "unstructured".
  measurement_error : bool, optional
      Whether or not to assume the endogenous observations `endog` were
      measured with error. Default is False.
  enforce_stationarity : bool, optional
      Whether or not to transform the AR parameters to enforce stationarity
      in the autoregressive component of the model. Default is True.
  enforce_invertibility : bool, optional
      Whether or not to transform the MA parameters to enforce invertibility
      in the moving average component of the model. Default is True.
  trend_offset : int, optional
      The offset at which to start time trend values. Default is 1, so that
      if `trend='t'` the trend is equal to 1, 2, ..., nobs. Typically is only
      set when the model created by extending a previous dataset.
  **kwargs
      Keyword arguments may be used to provide default values for state space
      matrices or for Kalman filtering options. See `Representation`, and
      `KalmanFilter` for more details.
  
  Attributes
  ----------
  order : iterable
      The (p,q) order of the model for the number of AR and MA parameters to
      use.
  trend : str{'n','c','t','ct'} or iterable
      Parameter controlling the deterministic trend polynomial :math:`A(t)`.
      Can be specified as a string where 'c' indicates a constant (i.e. a
      degree zero component of the trend polynomial), 't' indicates a
      linear trend with time, and 'ct' is both. Can also be specified as an
      iterable defining the non-zero polynomial exponents to include, in
      increasing order. For example, `[1,1,0,1]` denotes
      :math:`a + bt + ct^3`.
  error_cov_type : {'diagonal', 'unstructured'}, optional
      The structure of the covariance matrix of the error term, where
      "unstructured" puts no restrictions on the matrix and "diagonal"
      requires it to be a diagonal matrix (uncorrelated errors). Default is
      "unstructured".
  measurement_error : bool, optional
      Whether or not to assume the endogenous observations `endog` were
      measured with error. Default is False.
  enforce_stationarity : bool, optional
      Whether or not to transform the AR parameters to enforce stationarity
      in the autoregressive component of the model. Default is True.
  enforce_invertibility : bool, optional
      Whether or not to transform the MA parameters to enforce invertibility
      in the moving average component of the model. Default is True.
  
  Notes
  -----
  Generically, the VARMAX model is specified (see for example chapter 18 of
  [1]_):
  
  .. math::
  
      y_t = A(t) + A_1 y_{t-1} + \dots + A_p y_{t-p} + B x_t + \epsilon_t +
      M_1 \epsilon_{t-1} + \dots M_q \epsilon_{t-q}
  
  where :math:`\epsilon_t \sim N(0, \Omega)`, and where :math:`y_t` is a
  `k_endog x 1` vector. Additionally, this model allows considering the case
  where the variables are measured with error.
  
  Note that in the full VARMA(p,q) case there is a fundamental identification
  problem in that the coefficient matrices :math:`\{A_i, M_j\}` are not
  generally unique, meaning that for a given time series process there may
  be multiple sets of matrices that equivalently represent it. See Chapter 12
  of [1]_ for more information. Although this class can be used to estimate
  VARMA(p,q) models, a warning is issued to remind users that no steps have
  been taken to ensure identification in this case.
  
  References
  ----------
  .. [1] Lütkepohl, Helmut. 2007.
     New Introduction to Multiple Time Series Analysis.
     Berlin: Springer.
code: def __init__(
      self,
      endog,
      exog=None,
      order=(1, 0),
      trend="c",
      error_cov_type="unstructured",
      measurement_error=False,
      enforce_stationarity=True,
      enforce_invertibility=True,
      trend_offset=1,
      **kwargs,
  ):
      # Model parameters
      self.error_cov_type = error_cov_type
      self.measurement_error = measurement_error
      self.enforce_stationarity = enforce_stationarity
      self.enforce_invertibility = enforce_invertibility
  
      # Save the given orders
      self.order = order
  
      # Model orders
      self.k_ar = int(order[0])
      self.k_ma = int(order[1])
  
      # Check for valid model
      if error_cov_type not in ["diagonal", "unstructured"]:
          raise ValueError(
              "Invalid error covariance matrix type" " specification."
          )
      if self.k_ar == 0 and self.k_ma == 0:
          raise ValueError(
              "Invalid VARMAX(p,q) specification; at least one"
              " p,q must be greater than zero."
          )
  
      # Warn for VARMA model
      if self.k_ar > 0 and self.k_ma > 0:
          warn(
              "Estimation of VARMA(p,q) models is not generically robust,"
              " due especially to identification issues.",
              EstimationWarning,
          )
  
      # Trend
      self.trend = trend
      self.trend_offset = trend_offset
      self.polynomial_trend, self.k_trend = prepare_trend_spec(
          self.trend
      )
      self._trend_is_const = (
          self.polynomial_trend.size == 1
          and self.polynomial_trend[0] == 1
      )
  
      # Exogenous data
      (self.k_exog, exog) = prepare_exog(exog)
  
      # Note: at some point in the future might add state regression, as in
      # SARIMAX.
      self.mle_regression = self.k_exog > 0
  
      # We need to have an array or pandas at this point
      if not _is_using_pandas(endog, None):
          endog = np.asanyarray(endog)
  
      # Model order
      # Used internally in various places
      _min_k_ar = max(self.k_ar, 1)
      self._k_order = _min_k_ar + self.k_ma
  
      # Number of states
      k_endog = endog.shape[1]
      k_posdef = k_endog
      k_states = k_endog * self._k_order
  
      # By default, initialize as stationary
      kwargs.setdefault("initialization", "stationary")
  
      # By default, use LU decomposition
      kwargs.setdefault(
          "inversion_method", INVERT_UNIVARIATE | SOLVE_LU
      )
  
      # Initialize the state space model
      super().__init__(
          endog,
          exog=exog,
          k_states=k_states,
          k_posdef=k_posdef,
          **kwargs,
      )
  
      # Set as time-varying model if we have time-trend or exog
      if self.k_exog > 0 or (
          self.k_trend > 0 and not self._trend_is_const
      ):
          self.ssm._time_invariant = False
  
      # Initialize the parameters
      self.parameters = {}
      self.parameters["trend"] = self.k_endog * self.k_trend
      self.parameters["ar"] = self.k_endog**2 * self.k_ar
      self.parameters["ma"] = self.k_endog**2 * self.k_ma
      self.parameters["regression"] = self.k_endog * self.k_exog
      if self.error_cov_type == "diagonal":
          self.parameters["state_cov"] = self.k_endog
      # These parameters fill in a lower-triangular matrix which is then
      # dotted with itself to get a positive definite matrix.
      elif self.error_cov_type == "unstructured":
          self.parameters["state_cov"] = int(
              self.k_endog * (self.k_endog + 1) / 2
          )
      self.parameters["obs_cov"] = self.k_endog * self.measurement_error
      self.k_params = sum(self.parameters.values())
  
      # Initialize trend data: we create trend data with one more observation
      # than we actually have, to make it easier to insert the appropriate
      # trend component into the final state intercept.
      trend_data = prepare_trend_data(
          self.polynomial_trend,
          self.k_trend,
          self.nobs + 1,
          offset=self.trend_offset,
      )
      self._trend_data = trend_data[:-1]
      self._final_trend = trend_data[-1:]
  
      # Initialize known elements of the state space matrices
  
      # If we have exog effects, then the state intercept needs to be
      # time-varying
      if (
          self.k_trend > 0 and not self._trend_is_const
      ) or self.k_exog > 0:
          self.ssm["state_intercept"] = np.zeros(
              (self.k_states, self.nobs)
          )
          # self.ssm['obs_intercept'] = np.zeros((self.k_endog, self.nobs))
  
      # The design matrix is just an identity for the first k_endog states
      idx = np.diag_indices(self.k_endog)
      self.ssm[("design",) + idx] = 1
  
      # The transition matrix is described in four blocks, where the upper
      # left block is in companion form with the autoregressive coefficient
      # matrices (so it is shaped k_endog * k_ar x k_endog * k_ar) ...
      if self.k_ar > 0:
          idx = np.diag_indices((self.k_ar - 1) * self.k_endog)
          idx = idx[0] + self.k_endog, idx[1]
          self.ssm[("transition",) + idx] = 1
      # ... and the  lower right block is in companion form with zeros as the
      # coefficient matrices (it is shaped k_endog * k_ma x k_endog * k_ma).
      idx = np.diag_indices((self.k_ma - 1) * self.k_endog)
      idx = (
          idx[0] + (_min_k_ar + 1) * self.k_endog,
          idx[1] + _min_k_ar * self.k_endog,
      )
      self.ssm[("transition",) + idx] = 1
  
      # The selection matrix is described in two blocks, where the upper
      # block selects the all k_posdef errors in the first k_endog rows
      # (the upper block is shaped k_endog * k_ar x k) and the lower block
      # also selects all k_posdef errors in the first k_endog rows (the lower
      # block is shaped k_endog * k_ma x k).
      idx = np.diag_indices(self.k_endog)
      self.ssm[("selection",) + idx] = 1
      idx = idx[0] + _min_k_ar * self.k_endog, idx[1]
      if self.k_ma > 0:
          self.ssm[("selection",) + idx] = 1
  
      # Cache some indices
      if self._trend_is_const and self.k_exog == 0:
          self._idx_state_intercept = np.s_[
              "state_intercept", :k_endog, :
          ]
      elif self.k_trend > 0 or self.k_exog > 0:
          self._idx_state_intercept = np.s_[
              "state_intercept", :k_endog, :-1
          ]
      if self.k_ar > 0:
          self._idx_transition = np.s_["transition", :k_endog, :]
      else:
          self._idx_transition = np.s_["transition", :k_endog, k_endog:]
      if self.error_cov_type == "diagonal":
          self._idx_state_cov = ("state_cov",) + np.diag_indices(
              self.k_endog
          )
      elif self.error_cov_type == "unstructured":
          self._idx_lower_state_cov = np.tril_indices(self.k_endog)
      if self.measurement_error:
          self._idx_obs_cov = ("obs_cov",) + np.diag_indices(
              self.k_endog
          )
  
      # Cache some slices
      def _slice(key, offset):
          length = self.parameters[key]
          param_slice = np.s_[offset : offset + length]
          offset += length
          return param_slice, offset
  
      offset = 0
      self._params_trend, offset = _slice("trend", offset)
      self._params_ar, offset = _slice("ar", offset)
      self._params_ma, offset = _slice("ma", offset)
      self._params_regression, offset = _slice("regression", offset)
      self._params_state_cov, offset = _slice("state_cov", offset)
      self._params_obs_cov, offset = _slice("obs_cov", offset)
  
      # Variable holding optional final `exog`
      # (note: self._final_trend was set earlier)
      self._final_exog = None
  
      # Update _init_keys attached by super
      self._init_keys += [
          "order",
          "trend",
          "error_cov_type",
          "measurement_error",
          "enforce_stationarity",
          "enforce_invertibility",
      ] + list(kwargs.keys())
  
imports: import contextlib
  from warnings import warn
  import pandas
  import numpy
  from statsmodels.compat.pandas import Appender
  from statsmodels.tools.tools import Bunch
  from statsmodels.tools.data import _is_using_pandas
  from statsmodels.tsa.vector_ar import var_model
  import statsmodels.base.wrapper
  from statsmodels.tools.sm_exceptions import EstimationWarning
  from .kalman_filter import INVERT_UNIVARIATE, SOLVE_LU
  from .mlemodel import MLEModel, MLEResults, MLEResultsWrapper
  from .initialization import Initialization
  from .tools import is_invertible, concat, prepare_exog, constrain_stationary_multivariate, unconstrain_stationary_multivariate, prepare_trend_spec, prepare_trend_data
label_desc: Use image node_id 1 to create a new VARMAX object from inherited base classes: MLEModel with example: obj = VARMAX(endog, exog, order, trend, error_cov_type, measurement_error, enforce_stationarity, enforce_invertibility, trend_offset)