Skip to content

Perturbation

Compute first order approximation of optimal controls

Parameters:

NumericModel

Model to be solved

boolean

If True: displays number of contracting eigenvalues

ndarray

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

boolean

Use nonlinear solver to find the steady-state

{1}

Approximation order. (Currently, only first order is supported).

Returns:

TaylorExpansion

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

Source code in src/dolo/algos/perturbation.py
def perturb(
    model: Model,
    *,
    details: bool = True,  #
    verbose: bool = True,  #
    steady_state=None,
    eigmax=1.0 - 1e-6,
    solve_steady_state=False,
    order=1,
):
    """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, IIDProcess
    from dolo.numeric.distribution import MvNormal

    process = model.exogenous

    if isinstance(process, VAR1):
        C_m = C[:, : len(m)]
        C_s = C[:, len(m) :]
    elif isinstance(process, IIDProcess):
        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
        )