docs: def geometric_sums(self, beta, x_t):
          r"""
          Forecast the geometric sums
  
          .. math::
  
              S_x := E \Big[ \sum_{j=0}^{\infty} \beta^j x_{t+j} | x_t \Big]
  
              S_y := E \Big[ \sum_{j=0}^{\infty} \beta^j y_{t+j} | x_t \Big]
  
          Parameters
          ----------
          beta : scalar(float)
              Discount factor, in [0, 1)
  
          beta : array_like(float)
              The term x_t for conditioning
  
          Returns
          -------
          S_x : array_like(float)
              Geometric sum as defined above
  
          S_y : array_like(float)
              Geometric sum as defined above
  
          """
  
          I = np.identity(self.n)
          S_x = solve(I - beta * self.A, x_t)
          S_y = self.G.dot(S_x)
  
          return S_x, S_y