Compute first order approximation of optimal controls

    Parameters:
    -----------

    model: NumericModel
        Model to be solved

    verbose: boolean
        If True: displays number of contracting eigenvalues

    steady_state: ndarray
        Use supplied steady-state value to compute the approximation.
        The routine doesn't check whether it is really a solution or not.

    solve_steady_state: boolean
        Use nonlinear solver to find the steady-state

    orders: {1}
        Approximation order. (Currently, only first order is supported).

    Returns:
    --------

    TaylorExpansion:
        Decision Rule for the optimal controls around the steady-state.

def perturb(model, verbose=False, steady_state=None, eigmax=1.0-1e-6,
                         solve_steady_state=False, order=1, details=True):
    """Compute first order approximation of optimal controls

    Parameters:
    -----------

    model: NumericModel
        Model to be solved

    verbose: boolean
        If True: displays number of contracting eigenvalues

    steady_state: ndarray
        Use supplied steady-state value to compute the approximation.
        The routine doesn't check whether it is really a solution or not.

    solve_steady_state: boolean
        Use nonlinear solver to find the steady-state

    orders: {1}
        Approximation order. (Currently, only first order is supported).

    Returns:
    --------

    TaylorExpansion:
        Decision Rule for the optimal controls around the steady-state.

    """

    if order > 1:
        raise Exception("Not implemented.")

    if steady_state is None:
        steady_state = model.calibration

    G_s, G_x, G_e, F_s, F_x, F_S, F_X = get_derivatives(model, steady_state=steady_state)

    C, eigvals = approximate_1st_order(G_s, G_x, G_e, F_s, F_x, F_S, F_X)

    m = steady_state['exogenous']
    s = steady_state['states']
    x = steady_state['controls']

    from dolo.numeric.processes import VAR1
    from dolo.numeric.processes import MvNormal

    process = model.exogenous

    if isinstance(process, VAR1):
        C_m = C[:,:len(m)]
        C_s = C[:,len(m):]
    elif isinstance(process, MvNormal):
        C_m = None
        C_s = C
    dr = BivariateTaylor(m,s,x,C_m,C_s)
    if not details:
        return dr
    else:
        return PerturbationResult(
            dr,
            eigvals,
            True, # otherwise an Exception should have been raised already
            True, # otherwise an Exception should have been raised already
            True  # otherwise an Exception should have been raised already
        )