docs: def robust_rule(self, method='doubling'):
          """
          This method solves the robust control problem by tricking it
          into a stacked LQ problem, as described in chapter 2 of Hansen-
          Sargent's text "Robustness."  The optimal control with observed
          state is
  
          .. math::
  
              u_t = - F x_t
  
          And the value function is :math:`-x'Px`
  
          Parameters
          ----------
          method : str, optional(default='doubling')
              Solution method used in solving the associated Riccati
              equation, str in {'doubling', 'qz'}.
              
          Returns
          -------
          F : array_like(float, ndim=2)
              The optimal control matrix from above
          P : array_like(float, ndim=2)
              The positive semi-definite matrix defining the value
              function
          K : array_like(float, ndim=2)
              the worst-case shock matrix K, where
              :math:`w_{t+1} = K x_t` is the worst case shock
  
          """
          # == Simplify names == #
          A, B, C, Q, R = self.A, self.B, self.C, self.Q, self.R
          beta, theta = self.beta, self.theta
          k, j = self.k, self.j
          # == Set up LQ version == #
          I = identity(j)
          Z = np.zeros((k, j))
  
          if self.pure_forecasting:
              lq = LQ(-beta*I*theta, R, A, C, beta=beta)
  
              # == Solve and convert back to robust problem == #
              P, f, d = lq.stationary_values(method=method)
              F = np.zeros((self.k, self.n))
              K = -f[:k, :]
  
          else:
              Ba = hstack([B, C])
              Qa = vstack([hstack([Q, Z]), hstack([Z.T, -beta*I*theta])])
              lq = LQ(Qa, R, A, Ba, beta=beta)
  
              # == Solve and convert back to robust problem == #
              P, f, d = lq.stationary_values(method=method)
              F = f[:k, :]
              K = -f[k:f.shape[0], :]
  
          return F, K, P