How to use the quantecon.LQ function in quantecon

"""
    Taking the parameters A, B, Q, R as found in the `setup_matrices`,
    we find the value function of the optimal linear regulator problem.
    This is steps 2 and 3 in the lecture notes.
    Parameters
    ----------
    (A, B, Q, R) : Array(Float, ndim=2)
        The matrices that describe the oligopoly problem
    Returns
    -------
    (P, F, d) : Array(Float, ndim=2)
        The matrix that describes the value function of the optimal
        linear regulator problem.
    """

    lq = LQ(Q, -R, A, B, beta=beta)
    P, F, d = lq.stationary_values()

    Af = np.vstack((np.hstack([A-np.dot(B,F), np.array([[0., 0., 0., 0., 0.]]).T]),np.array([[0., 0., 0., 0., 0., 1.]])))
    Bf = np.array([[0., 0., 0., 0., 0., 1.]]).T

    lqf = LQ(Q, -Rf, Af, Bf, beta=beta)
    Pf, Ff, df = lqf.stationary_values()

    return P, F, d, Pf, Ff, df

R = np.array([[     0,  -A0 / 2,      0,      0],
                 [-A0 / 2,   A1 / 2, -μ / 2,      0],
                 [      0,   -μ / 2,      0,      0],
                 [      0,        0,      0,  d / 2]])

    A = np.array([[     1,      0,   0,              0],
                 [      0,      1,   0,              1],
                 [      0,      0,   0,              0],
                 [-A0 / d, A1 / d,   0, A1 / d + 1 / β]])

    B = np.array([0, 0, 1, 1/d]).reshape(-1, 1)

    Q = 0

    # Use LQ to solve the Ramsey Problem.
    lq = LQ(Q, -R, A, B, beta=β)
    P, F, d = lq.stationary_values()

    # Need y_0 to compute government tax revenue.
    P21 = P[3, :3]
    P22 = P[3, 3]
    z0 = np.array([1, Q0, τ0]).reshape(-1, 1)
    u0 = -P22**(-1) * P21 @ z0
    y0 = np.vstack([z0, u0])

    # Define A_F and S matricies
    AF = A - B @ F
    S = np.array([0, 1, 0, 0]).reshape(-1, 1) @ np.array([[0, 0, 1, 0]])

    # Solves equation (25)
    temp = β * AF.T @ S @ AF
    Ω = solve_discrete_lyapunov(np.sqrt(β) * AF.T, temp)

"""
    Taking the parameters A, B, Q, R as found in the `setup_matrices`,
    we find the value function of the optimal linear regulator problem.
    This is steps 2 and 3 in the lecture notes.
    Parameters
    ----------
    (A, B, Q, R) : Array(Float, ndim=2)
        The matrices that describe the oligopoly problem
    Returns
    -------
    (P, F, d) : Array(Float, ndim=2)
        The matrix that describes the value function of the optimal
        linear regulator problem.
    """

    lq = LQ(Q, -R, A, B, beta=beta)
    P, F, d = lq.stationary_values()

    Af = np.vstack((np.hstack([A-np.dot(B,F), np.array([[0., 0., 0., 0., 0.]]).T]),np.array([[0., 0., 0., 0., 0., 1.]])))
    Bf = np.array([[0., 0., 0., 0., 0., 1.]]).T

    lqf = LQ(Q, -Rf, Af, Bf, beta=beta)
    Pf, Ff, df = lqf.stationary_values()

    return P, F, d, Pf, Ff, df

R = np.array([[0, -A0/2, 0, 0],
                 [-A0/2, A1/2, -mu/2, 0],
                 [0, -mu/2, 0, 0],
                 [0, 0, 0, d/2]])

    A = np.array([[1, 0, 0, 0],
                 [0, 1, 0, 1],
                 [0, 0, 0, 0],
                 [-A0/d, A1/d, 0, A1/d+1/beta]])

    B = np.array([0, 0, 1, 1/d]).reshape(-1, 1)

    Q = 0

    # Use LQ to solve the Ramsey Problem.
    lq = LQ(Q, -R, A, B, beta=beta)
    P, F, d = lq.stationary_values()

    # Need y_0 to compute government tax revenue.
    P21 = P[3, :3]
    P22 = P[3, 3]
    z0 = np.array([1, Q0, tau0]).reshape(-1, 1)
    u0 = -P22**(-1) * P21.dot(z0)
    y0 = np.vstack([z0, u0])

    # Define A_F and S matricies
    AF = A - B.dot(F)
    S = np.array([0, 1, 0, 0]).reshape(-1, 1).dot(np.array([[0, 0, 1, 0]]))

    # Solves equation (25)
    temp = beta * AF.T.dot(S).dot(AF)
    Omega = solve_discrete_lyapunov(np.sqrt(beta) * AF.T, temp)

(A, B, Q, R) : Array(Float, ndim=2)
        The matrices that describe the oligopoly problem
    Returns
    -------
    (P, F, d) : Array(Float, ndim=2)
        The matrix that describes the value function of the optimal
        linear regulator problem.
    """

    lq = LQ(Q, -R, A, B, beta=beta)
    P, F, d = lq.stationary_values()

    Af = np.vstack((np.hstack([A-np.dot(B,F), np.array([[0., 0., 0., 0., 0.]]).T]),np.array([[0., 0., 0., 0., 0., 1.]])))
    Bf = np.array([[0., 0., 0., 0., 0., 1.]]).T

    lqf = LQ(Q, -Rf, Af, Bf, beta=beta)
    Pf, Ff, df = lqf.stationary_values()

    return P, F, d, Pf, Ff, df

(A, B, Q, R) : Array(Float, ndim=2)
        The matrices that describe the oligopoly problem
    Returns
    -------
    (P, F, d) : Array(Float, ndim=2)
        The matrix that describes the value function of the optimal
        linear regulator problem.
    """

    lq = LQ(Q, -R, A, B, beta=beta)
    P, F, d = lq.stationary_values()

    Af = np.vstack((np.hstack([A-np.dot(B,F), np.array([[0., 0., 0., 0., 0.]]).T]),np.array([[0., 0., 0., 0., 0., 1.]])))
    Bf = np.array([[0., 0., 0., 0., 0., 1.]]).T

    lqf = LQ(Q, -Rf, Af, Bf, beta=beta)
    Pf, Ff, df = lqf.stationary_values()

    return P, F, d, Pf, Ff, df

from quantecon import LQ
from solution_ree_ex1 import beta, R, Q, B

candidates = (
          (94.0886298678, 0.923409232937),
          (93.2119845412, 0.984323478873),
          (95.0818452486, 0.952459076301)
             )

for kappa0, kappa1 in candidates:

    # == Form the associated law of motion == #
    A = np.array([[1, 0, 0], [0, kappa1, kappa0], [0, 0, 1]])

    # == Solve the LQ problem for the firm == #
    lq = LQ(Q, R, A, B, beta=beta)
    P, F, d = lq.stationary_values()
    F = F.flatten()
    h0, h1, h2 = -F[2], 1 - F[0], -F[1]

    # == Test the equilibrium condition == #
    if np.allclose((kappa0, kappa1), (h0, h1 + h2)):
        print('Equilibrium pair =', kappa0, kappa1)
        print('(h0, h1, h2) = ', h0, h1, h2)
        break