docs: def stationary_values(self, method='doubling'):
          """
          Computes the matrix :math:`P` and scalar :math:`d` that represent
          the value function
  
          .. math::
  
               V(x) = x' P x + d
  
          in the infinite horizon case.  Also computes the control matrix
          :math:`F` from :math:`u = - Fx`. Computation is via the solution
          algorithm as specified by the `method` option (default to the
          doubling algorithm) (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
          -------
          P : array_like(float)
              P is part of the value function representation of
              :math:`V(x) = x'Px + d`
          F : array_like(float)
              F is the policy rule that determines the choice of control
              in each period.
          d : array_like(float)
              d is part of the value function representation of
              :math:`V(x) = x'Px + d`
  
          """
          # === simplify notation === #
          Q, R, A, B, N, C = self.Q, self.R, self.A, self.B, self.N, self.C
  
          # === solve Riccati equation, obtain P === #
          A0, B0 = np.sqrt(self.beta) * A, np.sqrt(self.beta) * B
          P = solve_discrete_riccati(A0, B0, R, Q, N, method=method)
  
          # == Compute F == #
          S1 = Q + self.beta * dot(B.T, dot(P, B))
          S2 = self.beta * dot(B.T, dot(P, A)) + N
          F = solve(S1, S2)
  
          # == Compute d == #
          d = self.beta * np.trace(dot(P, dot(C, C.T))) / (1 - self.beta)
  
          # == Bind states and return values == #
          self.P, self.F, self.d = P, F, d
  
          return P, F, d