Computes the residuals, the derivatives of the stacked-time system.
    :param model: an fga model
    :param guess: the guess for the simulated values. An `(n_s.n_x) x N` array,
                  where n_s is the number of states,
    n_x the number of controls, and `N` the length of the simulation.
    :param start: initial boundary condition (initial value of the states)
    :param final: final boundary condition (last value of the controls)
    :param shocks: values for the exogenous shocks
    :param diff: if True, the derivatives are computes
    :return: a list with two elements:
        - an `(n_s.n_x) x N` array with the residuals of the system
        - a `(n_s.n_x) x N x (n_s.n_x) x N` array representing the jacobian of
             the system

def det_residual(model,
                 guess,
                 start,
                 final,
                 shocks,
                 diff=True,
                 jactype='sparse'):
    '''
    Computes the residuals, the derivatives of the stacked-time system.
    :param model: an fga model
    :param guess: the guess for the simulated values. An `(n_s.n_x) x N` array,
                  where n_s is the number of states,
    n_x the number of controls, and `N` the length of the simulation.
    :param start: initial boundary condition (initial value of the states)
    :param final: final boundary condition (last value of the controls)
    :param shocks: values for the exogenous shocks
    :param diff: if True, the derivatives are computes
    :return: a list with two elements:
        - an `(n_s.n_x) x N` array with the residuals of the system
        - a `(n_s.n_x) x N x (n_s.n_x) x N` array representing the jacobian of
             the system
    '''

    # TODO: compute a sparse derivative and ensure the solvers can deal with it

    n_s = len(model.symbols['states'])
    n_x = len(model.symbols['controls'])

    # n_e = len(model.symbols['shocks'])
    N = guess.shape[0]

    p = model.calibration['parameters']

    f = model.functions['arbitrage']
    g = model.functions['transition']

    vec = guess[:-1, :]
    vec_f = guess[1:, :]

    s = vec[:, :n_s]
    x = vec[:, n_s:]
    S = vec_f[:, :n_s]
    X = vec_f[:, n_s:]

    m = shocks[:-1, :]
    M = shocks[1:, :]

    if diff:
        SS, SS_m, SS_s, SS_x, SS_M = g(m, s, x, M, p, diff=True)
        R, R_m, R_s, R_x, R_M, R_S, R_X = f(m, s, x, M, S, X, p, diff=True)
    else:
        SS = g(m, s, x, M, p)
        R = f(m, s, x, M, S, X, p)

    res_s = SS - S
    res_x = R

    res = np.zeros((N, n_s + n_x))

    res[1:, :n_s] = res_s
    res[:-1, n_s:] = res_x

    res[0, :n_s] = -(guess[0, :n_s] - start)
    res[-1, n_s:] = -(guess[-1, n_s:] - guess[-2, n_s:])

    if not diff:
        return res
    else:

        sparse_jac = False
        if not sparse_jac:

            # we compute the derivative matrix
            res_s_s = SS_s
            res_s_x = SS_x

            # next block is probably very inefficient
            jac = np.zeros((N, n_s + n_x, N, n_s + n_x))
            for i in range(N - 1):
                jac[i, n_s:, i, :n_s] = R_s[i, :, :]
                jac[i, n_s:, i, n_s:] = R_x[i, :, :]
                jac[i, n_s:, i + 1, :n_s] = R_S[i, :, :]
                jac[i, n_s:, i + 1, n_s:] = R_X[i, :, :]
                jac[i + 1, :n_s, i, :n_s] = SS_s[i, :, :]
                jac[i + 1, :n_s, i, n_s:] = SS_x[i, :, :]
                jac[i + 1, :n_s, i + 1, :n_s] = -np.eye(n_s)
                # jac[i,n_s:,i,:n_s] = R_s[i,:,:]
                # jac[i,n_s:,i,n_s:] = R_x[i,:,:]
                # jac[i+1,n_s:,i,:n_s] = R_S[i,:,:]
                # jac[i+1,n_s:,i,n_s:] = R_X[i,:,:]
                # jac[i,:n_s,i+1,:n_s] = SS_s[i,:,:]
                # jac[i,:n_s,i+1,n_s:] = SS_x[i,:,:]
                # jac[i+1,:n_s,i+1,:n_s] = -np.eye(n_s)
            jac[0, :n_s, 0, :n_s] = -np.eye(n_s)
            jac[-1, n_s:, -1, n_s:] = -np.eye(n_x)
            jac[-1, n_s:, -2, n_s:] = +np.eye(n_x)
            nn = jac.shape[0] * jac.shape[1]
            res = res.ravel()
            jac = jac.reshape((nn, nn))

        if jactype == 'sparse':
            from scipy.sparse import csc_matrix, csr_matrix
            jac = csc_matrix(jac)
            # scipy bug ? I don't get the same with csr

        return [res, jac]