docs: def compute_sequence(self, x0, ts_length=None, method='doubling',
                           random_state=None):
          """
          Compute and return the optimal state and control sequences
          :math:`x_0, ..., x_T` and :math:`u_0,..., u_T`  under the
          assumption that :math:`{w_t}` is iid and :math:`N(0, 1)`.
  
          Parameters
          ----------
          x0 : array_like(float)
              The initial state, a vector of length n
  
          ts_length : scalar(int)
              Length of the simulation -- defaults to T in finite case
  
          method : str, optional(default='doubling')
              Solution method used in solving the associated Riccati
              equation, str in {'doubling', 'qz'}. Only relevant when the
              `T` attribute is `None` (i.e., the horizon is infinite).
  
          random_state : int or np.random.RandomState, optional
              Random seed (integer) or np.random.RandomState instance to set
              the initial state of the random number generator for
              reproducibility. If None, a randomly initialized RandomState is
              used.
  
          Returns
          -------
          x_path : array_like(float)
              An n x T+1 matrix, where the t-th column represents :math:`x_t`
  
          u_path : array_like(float)
              A k x T matrix, where the t-th column represents :math:`u_t`
  
          w_path : array_like(float)
              A j x T+1 matrix, where the t-th column represent :math:`w_t`
  
          """
  
          # === Simplify notation === #
          A, B, C = self.A, self.B, self.C
  
          # == Preliminaries, finite horizon case == #
          if self.T:
              T = self.T if not ts_length else min(ts_length, self.T)
              self.P, self.d = self.Rf, 0
  
          # == Preliminaries, infinite horizon case == #
          else:
              T = ts_length if ts_length else 100
              self.stationary_values(method=method)
  
          # == Set up initial condition and arrays to store paths == #
          random_state = check_random_state(random_state)
          x0 = np.asarray(x0)
          x0 = x0.reshape(self.n, 1)  # Make sure x0 is a column vector
          x_path = np.empty((self.n, T+1))
          u_path = np.empty((self.k, T))
          w_path = random_state.randn(self.j, T+1)
          Cw_path = dot(C, w_path)
  
          # == Compute and record the sequence of policies == #
          policies = []
          for t in range(T):
              if self.T:  # Finite horizon case
                  self.update_values()
              policies.append(self.F)
  
          # == Use policy sequence to generate states and controls == #
          F = policies.pop()
          x_path[:, 0] = x0.flatten()
          u_path[:, 0] = - dot(F, x0).flatten()
          for t in range(1, T):
              F = policies.pop()
              Ax, Bu = dot(A, x_path[:, t-1]), dot(B, u_path[:, t-1])
              x_path[:, t] = Ax + Bu + Cw_path[:, t]
              u_path[:, t] = - dot(F, x_path[:, t])
          Ax, Bu = dot(A, x_path[:, T-1]), dot(B, u_path[:, T-1])
          x_path[:, T] = Ax + Bu + Cw_path[:, T]
  
          return x_path, u_path, w_path