docs: def compute_sequence(self, x0, ts_length=None, Pay=None):
          """
          Simulate quantities and prices for the economy
  
          Parameters
          ----------
          x0 : array_like(float)
              The initial state
  
          ts_length : scalar(int)
              Length of the simulation
  
          Pay : array_like(float)
              Vector to price an asset whose payout is Pay*xt
  
          """
          lq = LQ(self.Q, self.R, self.A, self.B,
                  self.C, N=self.W, beta=self.beta)
          xp, up, wp = lq.compute_sequence(x0, ts_length)
          self.h = self.Sh.dot(xp)
          self.k = self.Sk.dot(xp)
          self.i = self.Si.dot(xp)
          self.b = self.Sb.dot(xp)
          self.d = self.Sd.dot(xp)
          self.c = self.Sc.dot(xp)
          self.g = self.Sg.dot(xp)
          self.s = self.Ss.dot(xp)
  
          # === Value of J-period risk-free bonds === #
          # === See p.145: Equation (7.11.2) === #
          e1 = np.zeros((1, self.nc))
          e1[0, 0] = 1
          self.R1_Price = np.empty((ts_length + 1, 1))
          self.R2_Price = np.empty((ts_length + 1, 1))
          self.R5_Price = np.empty((ts_length + 1, 1))
          for i in range(ts_length + 1):
              self.R1_Price[i, 0] = self.beta * e1.dot(self.Mc).dot(np.linalg.matrix_power(
                  self.A0, 1)).dot(xp[:, i]) / e1.dot(self.Mc).dot(xp[:, i])
              self.R2_Price[i, 0] = self.beta**2 * e1.dot(self.Mc).dot(
                  np.linalg.matrix_power(self.A0, 2)).dot(xp[:, i]) / e1.dot(self.Mc).dot(xp[:, i])
              self.R5_Price[i, 0] = self.beta**5 * e1.dot(self.Mc).dot(
                  np.linalg.matrix_power(self.A0, 5)).dot(xp[:, i]) / e1.dot(self.Mc).dot(xp[:, i])
  
          # === Gross rates of return on 1-period risk-free bonds === #
          self.R1_Gross = 1 / self.R1_Price
  
          # === Net rates of return on J-period risk-free bonds === #
          # === See p.148: log of gross rate of return, divided by j === #
          self.R1_Net = np.log(1 / self.R1_Price) / 1
          self.R2_Net = np.log(1 / self.R2_Price) / 2
          self.R5_Net = np.log(1 / self.R5_Price) / 5
  
          # === Value of asset whose payout vector is Pay*xt === #
          # See p.145: Equation (7.11.1)
          if isinstance(Pay, np.ndarray) == True:
              self.Za = Pay.T.dot(self.Mc)
              self.Q = solve_discrete_lyapunov(
                  self.A0.T * self.beta**0.5, self.Za)
              self.q = self.beta / (1 - self.beta) * \
                  np.trace(self.C.T.dot(self.Q).dot(self.C))
              self.Pay_Price = np.empty((ts_length + 1, 1))
              self.Pay_Gross = np.empty((ts_length + 1, 1))
              self.Pay_Gross[0, 0] = np.nan
              for i in range(ts_length + 1):
                  self.Pay_Price[i, 0] = (xp[:, i].T.dot(self.Q).dot(
                      xp[:, i]) + self.q) / e1.dot(self.Mc).dot(xp[:, i])
              for i in range(ts_length):
                  self.Pay_Gross[i + 1, 0] = self.Pay_Price[i + 1,
                                                            0] / (self.Pay_Price[i, 0] - Pay.dot(xp[:, i]))
          return