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

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