content: """
  Statespace Tools
  
  Author: Chad Fulton
  License: Simplified-BSD
  """
  from __future__ import division, absolute_import, print_function
  
  import numpy as np
  from statsmodels.tools.data import _is_using_pandas
  from . import _statespace
  
  try:
      from scipy.linalg.blas import find_best_blas_type
  except ImportError:  # pragma: no cover
      # Shim for SciPy 0.11, derived from tag=0.11 scipy.linalg.blas
      _type_conv = {'f': 's', 'd': 'd', 'F': 'c', 'D': 'z', 'G': 'z'}
  
      def find_best_blas_type(arrays):
          dtype, index = max(
              [(ar.dtype, i) for i, ar in enumerate(arrays)])
          prefix = _type_conv.get(dtype.char, 'd')
          return prefix, dtype, None
  
  
  prefix_dtype_map = {
      's': np.float32, 'd': np.float64, 'c': np.complex64, 'z': np.complex128
  }
  prefix_statespace_map = {
      's': _statespace.sStatespace, 'd': _statespace.dStatespace,
      'c': _statespace.cStatespace, 'z': _statespace.zStatespace
  }
  prefix_kalman_filter_map = {
      's': _statespace.sKalmanFilter, 'd': _statespace.dKalmanFilter,
      'c': _statespace.cKalmanFilter, 'z': _statespace.zKalmanFilter
  }
  
  
  def companion_matrix(polynomial):
      r"""
      Create a companion matrix
  
      Parameters
      ----------
      polynomial : array_like, optional.
          If an iterable, interpreted as the coefficients of the polynomial from
          which to form the companion matrix. Polynomial coefficients are in
          order of increasing degree. If an integer, the size of the companion
          matrix (the polynomial coefficients are then set to zeros).
  
      Returns
      -------
      companion_matrix : array
  
      Notes
      -----
  
      Returns a matrix of the form
  
      .. math::
          \begin{bmatrix}
              \phi_1 & 1      & 0 & \cdots & 0 \\
              \phi_2 & 0      & 1 &        & 0 \\
              \vdots &        &   & \ddots & 0 \\
                     &        &   &        & 1 \\
              \phi_n & 0      & 0 & \cdots & 0 \\
          \end{bmatrix}
  
      where some or all of the :math:`\phi_i` may be non-zero (if `polynomial` is
      None, then all are equal to zero).
  
      If the coefficients provided are :math:`(c_0, c_1, \dots, c_{n})`,
      then the companion matrix is an :math:`n \times n` matrix formed with the
      elements in the first column defined as
      :math:`\phi_i = -\frac{c_i}{c_0}, i \in 1, \dots, n`.
      """
      if isinstance(polynomial, int):
          n = polynomial
          polynomial = None
      else:
          n = len(polynomial) - 1
          polynomial = np.asanyarray(polynomial)
  
      matrix = np.zeros((n, n))
      idx = np.diag_indices(n - 1)
      idx = (idx[0], idx[1] + 1)
      matrix[idx] = 1
      if polynomial is not None and n > 0:
          matrix[:, 0] = -polynomial[1:] / polynomial[0]
      return matrix
  
  
  def diff(series, k_diff=1, k_seasonal_diff=None, k_seasons=1):
      r"""
      Difference a series simply and/or seasonally along the zero-th axis.
  
      Given a series (denoted :math:`y_t`), performs the differencing operation
  
      .. math::
  
          \Delta^d \Delta_s^D y_t
  
      where :math:`d =` `diff`, :math:`s =` `k_seasons`,
      :math:`D =` `seasonal\_diff`, and :math:`\Delta` is the difference
      operator.
  
      Parameters
      ----------
      series : array_like
          The series to be differenced.
      diff : int, optional
          The number of simple differences to perform. Default is 1.
      seasonal_diff : int or None, optional
          The number of seasonal differences to perform. Default is no seasonal
          differencing.
      k_seasons : int, optional
          The seasonal lag. Default is 1. Unused if there is no seasonal
          differencing.
  
      Returns
      -------
      differenced : array
          The differenced array.
      """
      pandas = _is_using_pandas(series, None)
      differenced = np.asanyarray(series) if not pandas else series
  
      # Seasonal differencing
      if k_seasonal_diff is not None:
          while k_seasonal_diff > 0:
              if not pandas:
                  differenced = (
                      differenced[k_seasons:] - differenced[:-k_seasons]
                  )
              else:
                  differenced = differenced.diff(k_seasons)[k_seasons:]
              k_seasonal_diff -= 1
  
      # Simple differencing
      if not pandas:
          differenced = np.diff(differenced, k_diff, axis=0)
      else:
          while k_diff > 0:
              differenced = differenced.diff()[1:]
              k_diff -= 1
      return differenced
  
  
  def is_invertible(polynomial, threshold=1.):
      r"""
      Determine if a polynomial is invertible.
  
      Requires all roots of the polynomial lie inside the unit circle.
  
      Parameters
      ----------
      polynomial : array_like
          Coefficients of a polynomial, in order of increasing degree.
          For example, `polynomial=[1, -0.5]` corresponds to the polynomial
          :math:`1 - 0.5x` which has root :math:`2`.
      threshold : number
          Allowed threshold for `is_invertible` to return True. Default is 1.
  
      Notes
      -----
  
      If the coefficients provided are :math:`(c_0, c_1, \dots, c_n)`, then
      the corresponding polynomial is :math:`c_0 + c_1 L + \dots + c_n L^n`.
  
      There are three equivalent methods of determining if the polynomial
      represented by the coefficients is invertible:
  
      The first method factorizes the polynomial into:
  
      .. math::
  
          C(L) & = c_0 + c_1 L + \dots + c_n L^n \\
               & = constant (1 - \lambda_1 L)
                   (1 - \lambda_2 L) \dots (1 - \lambda_n L)
  
      In order for :math:`C(L)` to be invertible, it must be that each factor
      :math:`(1 - \lambda_i L)` is invertible; the condition is then that
      :math:`|\lambda_i| < 1`, where :math:`\lambda_i` is a root of the
      polynomial.
  
      The second method factorizes the polynomial into:
  
      .. math::
  
          C(L) & = c_0 + c_1 L + \dots + c_n L^n \\
               & = constant (L - \zeta_1 L) (L - \zeta_2) \dots (L - \zeta_3)
  
      The condition is now :math:`|\zeta_i| > 1`, where :math:`\zeta_i` is a root
      of the polynomial with reversed coefficients and
      :math:`\lambda_i = \frac{1}{\zeta_i}`.
  
      Finally, a companion matrix can be formed using the coefficients of the
      polynomial. Then the eigenvalues of that matrix give the roots of the
      polynomial. This last method is the one actually used.
  
      See Also
      --------
      companion_matrix
      """
      # First method:
      # np.all(np.abs(np.roots(np.r_[1, params])) < 1)
      # Second method:
      # np.all(np.abs(np.roots(np.r_[1, params][::-1])) > 1)
      # Final method:
      eigvals = np.linalg.eigvals(companion_matrix(polynomial))
      return np.all(np.abs(eigvals) < threshold)
  
  
  def constrain_stationary_univariate(unconstrained):
      """
      Transform unconstrained parameters used by the optimizer to constrained
      parameters used in likelihood evaluation
  
      Parameters
      ----------
      unconstrained : array
          Unconstrained parameters used by the optimizer, to be transformed to
          stationary coefficients of, e.g., an autoregressive or moving average
          component.
  
      Returns
      -------
      constrained : array
          Constrained parameters of, e.g., an autoregressive or moving average
          component, to be transformed to arbitrary parameters used by the
          optimizer.
  
      References
      ----------
  
      Monahan, John F. 1984.
      "A Note on Enforcing Stationarity in
      Autoregressive-moving Average Models."
      Biometrika 71 (2) (August 1): 403-404.
      """
  
      n = unconstrained.shape[0]
      y = np.zeros((n, n), dtype=unconstrained.dtype)
      r = unconstrained/((1 + unconstrained**2)**0.5)
      for k in range(n):
          for i in range(k):
              y[k, i] = y[k - 1, i] + r[k] * y[k - 1, k - i - 1]
          y[k, k] = r[k]
      return -y[n - 1, :]
  
  
  def unconstrain_stationary_univariate(constrained):
      """
      Transform constrained parameters used in likelihood evaluation
      to unconstrained parameters used by the optimizer
  
      Parameters
      ----------
      constrained : array
          Constrained parameters of, e.g., an autoregressive or moving average
          component, to be transformed to arbitrary parameters used by the
          optimizer.
  
      Returns
      -------
      unconstrained : array
          Unconstrained parameters used by the optimizer, to be transformed to
          stationary coefficients of, e.g., an autoregressive or moving average
          component.
  
      References
      ----------
  
      Monahan, John F. 1984.
      "A Note on Enforcing Stationarity in
      Autoregressive-moving Average Models."
      Biometrika 71 (2) (August 1): 403-404.
      """
      n = constrained.shape[0]
      y = np.zeros((n, n), dtype=constrained.dtype)
      y[n-1:] = -constrained
      for k in range(n-1, 0, -1):
          for i in range(k):
              y[k-1, i] = (y[k, i] - y[k, k]*y[k, k-i-1]) / (1 - y[k, k]**2)
      r = y.diagonal()
      x = r / ((1 - r**2)**0.5)
      return x
  
  
  def validate_matrix_shape(name, shape, nrows, ncols, nobs):
      """
      Validate the shape of a possibly time-varying matrix, or raise an exception
  
      Parameters
      ----------
      name : str
          The name of the matrix being validated (used in exception messages)
      shape : array_like
          The shape of the matrix to be validated. May be of size 2 or (if
          the matrix is time-varying) 3.
      nrows : int
          The expected number of rows.
      ncols : int
          The expected number of columns.
      nobs : int
          The number of observations (used to validate the last dimension of a
          time-varying matrix)
  
      Raises
      ------
      ValueError
          If the matrix is not of the desired shape.
      """
      ndim = len(shape)
  
      # Enforce dimension
      if ndim not in [2, 3]:
          raise ValueError('Invalid value for %s matrix. Requires a'
                           ' 2- or 3-dimensional array, got %d dimensions' %
                           (name, ndim))
      # Enforce the shape of the matrix
      if not shape[0] == nrows:
          raise ValueError('Invalid dimensions for %s matrix: requires %d'
                           ' rows, got %d' % (name, nrows, shape[0]))
      if not shape[1] == ncols:
          raise ValueError('Invalid dimensions for %s matrix: requires %d'
                           ' columns, got %d' % (name, ncols, shape[1]))
  
      # If we don't yet know `nobs`, don't allow time-varying arrays
      if nobs is None and not (ndim == 2 or shape[-1] == 1):
          raise ValueError('Invalid dimensions for %s matrix: time-varying'
                           ' matrices cannot be given unless `nobs` is specified'
                           ' (implicitly when a dataset is bound or else set'
                           ' explicity)' % name)
  
      # Enforce time-varying array size
      if ndim == 3 and nobs is not None and not shape[-1] in [1, nobs]:
          raise ValueError('Invalid dimensions for time-varying %s'
                           ' matrix. Requires shape (*,*,%d), got %s' %
                           (name, nobs, str(shape)))
  
  
  def validate_vector_shape(name, shape, nrows, nobs):
      """
      Validate the shape of a possibly time-varying vector, or raise an exception
  
      Parameters
      ----------
      name : str
          The name of the vector being validated (used in exception messages)
      shape : array_like
          The shape of the vector to be validated. May be of size 1 or (if
          the vector is time-varying) 2.
      nrows : int
          The expected number of rows (elements of the vector).
      nobs : int
          The number of observations (used to validate the last dimension of a
          time-varying vector)
  
      Raises
      ------
      ValueError
          If the vector is not of the desired shape.
      """
      ndim = len(shape)
      # Enforce dimension
      if ndim not in [1, 2]:
          raise ValueError('Invalid value for %s vector. Requires a'
                           ' 1- or 2-dimensional array, got %d dimensions' %
                           (name, ndim))
      # Enforce the shape of the vector
      if not shape[0] == nrows:
          raise ValueError('Invalid dimensions for %s vector: requires %d'
                           ' rows, got %d' % (name, nrows, shape[0]))
  
      # If we don't yet know `nobs`, don't allow time-varying arrays
      if nobs is None and not (ndim == 1 or shape[-1] == 1):
          raise ValueError('Invalid dimensions for %s vector: time-varying'
                           ' vectors cannot be given unless `nobs` is specified'
                           ' (implicitly when a dataset is bound or else set'
                           ' explicity)' % name)
  
      # Enforce time-varying array size
      if ndim == 2 and not shape[1] in [1, nobs]:
          raise ValueError('Invalid dimensions for time-varying %s'
                           ' vector. Requires shape (*,%d), got %s' %
                           (name, nobs, str(shape)))
  