Computes the limit of :math:`\Sigma_t` as t goes to infinity by
        solving the associated Riccati equation. The outputs are stored in the
        attributes `K_infinity` and `Sigma_infinity`. Computation is via the
        doubling algorithm (default) or a QZ decomposition method (see the
        documentation in `matrix_eqn.solve_discrete_riccati`).

        Parameters
        ----------
        method : str, optional(default="doubling")
            Solution method used in solving the associated Riccati
            equation, str in {'doubling', 'qz'}.

        Returns
        -------
        Sigma_infinity : array_like or scalar(float)
            The infinite limit of :math:`\Sigma_t`
        K_infinity : array_like or scalar(float)
            The stationary Kalman gain.

def stationary_values(self, method='doubling'):
        """
        Computes the limit of :math:`\Sigma_t` as t goes to infinity by
        solving the associated Riccati equation. The outputs are stored in the
        attributes `K_infinity` and `Sigma_infinity`. Computation is via the
        doubling algorithm (default) or a QZ decomposition method (see the
        documentation in `matrix_eqn.solve_discrete_riccati`).

        Parameters
        ----------
        method : str, optional(default="doubling")
            Solution method used in solving the associated Riccati
            equation, str in {'doubling', 'qz'}.

        Returns
        -------
        Sigma_infinity : array_like or scalar(float)
            The infinite limit of :math:`\Sigma_t`
        K_infinity : array_like or scalar(float)
            The stationary Kalman gain.

        """

        # === simplify notation === #
        A, C, G, H = self.ss.A, self.ss.C, self.ss.G, self.ss.H
        Q, R = np.dot(C, C.T), np.dot(H, H.T)

        # === solve Riccati equation, obtain Kalman gain === #
        Sigma_infinity = solve_discrete_riccati(A.T, G.T, Q, R, method=method)
        temp1 = dot(dot(A, Sigma_infinity), G.T)
        temp2 = inv(dot(G, dot(Sigma_infinity, G.T)) + R)
        K_infinity = dot(temp1, temp2)

        # == record as attributes and return == #
        self._Sigma_infinity, self._K_infinity = Sigma_infinity, K_infinity
        return Sigma_infinity, K_infinity