API Reference: Solvers & Solutions¶

def solve(
    A: TMatrix,
    B: TMatrix,
    C: TMatrix,
    method: Solver = "qz",
    options: dict | None = None,
) -> tuple[TMatrix, TVector | None]:
    """Solves AX² + BX + C = 0 for X using the chosen method

    Parameters
    ----------
    A, B, C : (N,N) Matrix

    method : str, optional
        chosen solver: either "ti" for fixed-point iteration or "qz" for generalized Schur decomposition, by default "qz"

    options : dict, optional
        dictionary of optional parameters to pass to the chosen solver, by default {}

    Returns
    -------
    (X, evs) : tuple[(N,N) Matrix, 2N Vector|None]
        solution of the equation as well as sorted list of associated generalized eigenvalues if the chosen method is "qz" and None otherwise

    Raises
    ------
    NoConvergence :
        when the convergence threshold is not reached within the maximal allowed number of iterations (only in `solve_ti`)
    LinAlgError :
        when a singular matrix is obtained during iterations in `solve_ti` or when Blachard-Kahn conditions are not verified in `solve_qz`
    ValueError :
        when a matrix containing a NaN is obtained
    """

    options = options or {}

    if method == "ti":
        sol, evs = solve_ti(A, B, C, **options)
    elif method == "qz":
        sol, evs = solve_qz(A, B, C, **options)
    else:
        raise ValueError(
            f"Unknown solver method: {method!r}. Valid choices are 'ti' and 'qz'."
        )

    return sol, evs

def solve_qz(
    A: TMatrix, B: TMatrix, C: TMatrix, tol: float = 1e-15
) -> tuple[TMatrix, TVector]:
    """Solves AX² + BX + C = 0 for X using QZ decomposition.

    Parameters
    ----------
    A, B, C : (N,N) Matrix

    tol : float, optional
        error tolerance, by default 1e-15

    Returns
    -------
    (X, evs) : tuple[(N,N) Matrix, 2N Vector]
        solution of the equation as well as sorted list of associated generalized eigenvalues

    Raises
    ------
    LinAlgError :
        when Blanchard-Kahn conditions are not verified (less than N generalized eigenvalues inside the unit ball)
    """
    from scipy.linalg import ordqz

    n = A.shape[0]
    if n == 0:
        return np.zeros((0, 0)), np.array([], dtype=float)
    I = np.eye(n)
    Z = np.zeros((n, n))

    # Generalised eigenvalue problem
    F = np.block([[Z, I], [-C, -B]])
    G = np.block([[I, Z], [Z, A]])
    tol_sort = 1e-6
    T, S, α, β, Q, Z = ordqz(F, G, sort=lambda a, b: np.abs(vgenev(a, b, tol=tol_sort)) <= 1 + tol_sort)  # type: ignore
    λ_all = vgenev(α, β, tol=tol_sort)
    Z11, Z12, Z21, Z22 = decompose_blocks(Z)

    λ_all = np.abs(λ_all)

    # TODO: verify whether Blanchard-Kahn conditions are valid
    evs = np.sort(λ_all).reshape(2 * n)
    n = len(evs) // 2
    l1 = evs[n - 1]
    l2 = evs[n]
    if l1 <= l2 < 1:
        raise BlanchardKahnError(
            f"Eigenvalue condition not satisfied: l_(n)={l1}, l_(n+1)={l2}. Too many stable solutions.",
            evs=evs,
        )
    if 1 < l1 <= l2:
        raise BlanchardKahnError(
            f"Eigenvalue condition not satisfied: l_(n)={l1}, l_(n+1)={l2}. No stable solution.",
            evs=evs,
        )

    X = (Z21 @ np.linalg.inv(Z11)).reshape(
        (n, n)
    )  # Reshape necessary for static type checking

    return X, evs

def solve_ti(
    A: TMatrix, B: TMatrix, C: TMatrix, T: int = 10000, tol: float = 1e-10
) -> tuple[TMatrix, None]:
    """Solves AX² + BX + C = 0 for X using fixed-point iteration.

    Parameters
    ----------
    A, B, C : (N,N) Matrix

    T : int, optional
        Maximum number of iterations. If more are needed, `NoConvergence` is raised, by default 10000

    tol : float, optional
        convergence threshold, by default 1e-10

    Returns
    -------
    (X, evs) : tuple[(N,N) Matrix, None]
        solution of the equation and None (necessary to have a common solver interface)

    Raises
    ------
    NoConvergence :
        when the convergence threshold is not reached within the maximal allowed number of iterations
    LinAlgError :
        when a singular matrix is obtained while iterating
    ValueError :
        when a matrix containing a NaN is obtained while iterating
    """
    n = A.shape[0]

    # Deterministic initial guess: the solution of the backward-looking part
    # (A = 0). Falls back to zeros when B is singular.
    X0: TMatrix
    try:
        X0 = np.asarray(linsolve(B, -C), dtype=float).reshape((n, n))
    except np.linalg.LinAlgError:
        X0 = np.zeros((n, n))

    for t in range(T):

        X1 = linsolve(A @ X0 + B, -C)
        e = abs(X0 - X1).max()

        if np.isnan(e):
            # impossible situation?
            raise ValueError("Invalid value")

        X0 = X1
        if e < tol:
            return X0, None

    raise NoConvergence("The maximal number of iterations was exceeded.")

def deterministic_solve(
    model,
    x0=None,
    T=None,
    method="hybr",
    verbose=False,
    continuation="stationary",
    growth_rate=None,
    growth_type="geometric",
    return_iterations=False,
    units: UnitsType = "level",
    **args,
):
    from .simul import TransitionSimulation

    continuation = args.pop("terminal_condition", continuation)
    growth_rate = args.pop("growth_rate", growth_rate)
    growth_type = args.pop("growth_type", growth_type)

    if x0 is None:
        v0 = model.deterministic_guess(T=T)
    else:
        v0 = np.array(x0)

    T = v0.shape[0] - 1

    u0 = np.array(v0).ravel()

    newton_res = newton(
        lambda u: model.deterministic_residuals_with_jacobian(
            u,
            sparsify=True,
            continuation=continuation,
            growth_rate=growth_rate,
            growth_type=growth_type,
            **args,
        ),
        u0,
        jactype="sparse",
        verbose=verbose,
        maxit=args.get("maxit", 10),
        tol=args.get("tol", 1e-8),
    )

    u_sol, nit = newton_res[0], newton_res[1]
    converged = getattr(newton_res, "converged", True)
    final_error = getattr(newton_res, "error", None)

    w0 = u_sol.reshape(v0.shape)

    variables = list(model.symbols["variables"])
    c = getattr(model, "context", {})
    ss_map = c.get("steady_states", {})
    val_map = c.get("values", {})
    ss_vec = np.array(
        [ss_map.get(name, val_map.get(name, {}).get(0, np.nan)) for name in variables],
        dtype=float,
    )

    attrs: dict[str, Any] = {
        "iterations": nit,
        "converged": converged,
    }
    if final_error is not None:
        attrs["residual"] = final_error

    sim = TransitionSimulation(
        data=w0,
        variables=variables,
        steady_state=ss_vec,
        units=units,
        canonical_units="level",
        attrs=attrs,
    )

    if return_iterations:
        return sim, nit

    return sim

def moments(X: TMatrix, Y: TMatrix, Σ: TMatrix) -> tuple[TMatrix, TMatrix]:
    """
    Computes conditional and unconditional moments of stationary process $y_t = X y_{t-1} + Y e_t$

    Parameters
    ----------
    X : (N,N) Matrix
        Transition matrix defining the stochastic process.
    Y : (N,N) Matrix
        Shock impact matrix defining the stochastic process.
    Sigma : (N,N) Matrix
        Covariance matrix of the independent identically distributed error terms e_t.

    Returns
    -------
    Gamma_0, Gamma : (N,N) Matrix
        Conditional and unconditional covariance matrices of the stationary process y_t respectively.

    Notes
    -----
    The unconditional covariance matrix Γ is computed in the following way:

    Applying the linear covariance operator to both sides of the equation $y_t = X y_{t-1} + Y e_t$ yields
    $$
    \\mathrm{Cov}(y_t) = X ⋅ \\mathrm{Cov}(y_{t-1}) ⋅ X^* + Y ⋅ \\mathrm{Cov}(e_t) ⋅ Y^*
    $$
    By stationarity of $y_t$, $\\mathrm{Cov}(y_t) = \\mathrm{Cov}(y_{t-1}) := Γ$, so
    $$
    Γ = X Γ X^* + Γ₀
    $$
    By applying the [Vec-operator](https://en.wikipedia.org/wiki/Vectorization_(mathematics)#Compatibility_with_Kronecker_products), we get the following equation:
    $$
    \\mathrm{Vec}(Γ) = (X ⊗ X) \\mathrm{Vec}(Γ)  + \\mathrm{Vec}(Γ₀)
    $$
    Which gives the following solution
    $$
    \\mathrm{Vec}(Γ) = (I_{N^2} - X ⊗ X)^{-1} \\mathrm{Vec}(Γ₀)
    $$
    """

    Γ0 = Y @ Σ @ Y.T
    n = X.shape[0]

    # Compute the unconditional variance
    Γ = (np.linalg.inv(np.eye(n**2) - np.kron(X, X)) @ Γ0.flatten()).reshape(n, n)

    return Γ0, Γ