'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
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`'

def compute_sequence(self, x0, ts_length=None):

    (A, B, C) = (self.A, self.B, self.C)
    if self.T:
        T = (self.T if (not ts_length) else min(ts_length, self.T))
        (self.P, self.d) = (self.Rf, 0)
    else:
        T = (ts_length if ts_length else 100)
        self.stationary_values()
    x0 = np.asarray(x0)
    x0 = x0.reshape(self.n, 1)
    x_path = np.empty((self.n, (T + 1)))
    u_path = np.empty((self.k, T))
    w_path = np.random.randn(self.j, (T + 1))
    Cw_path = dot(C, w_path)
    policies = []
    for t in range(T):
        if self.T:
            self.update_values()
        policies.append(self.F)
    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)
