Note:
a. Let $z_{t}=\left[\begin{array}{c} 1\\ s_{t}\\ w_{t} \end{array}\right]$ and $\epsilon_{t}=\left[\begin{array}{c} \epsilon_{1,t}\\ \epsilon_{2,t} \end{array}\right]$. Then, we can obtain
$$\begin{eqnarray} z_{t+1} &=& A_{22}z_{t}+C_{2}\epsilon_{t} \\ \left[\begin{array}{c} s_{t}\\ w_{t} \end{array}\right] &=& Sz_{t} \end{eqnarray}$$by letting $A_{22}=\left[\begin{array}{ccc} 1 & 0 & 0\\ \alpha_{s} & \rho_{1} & 0\\ \alpha_{w} & 0 & \gamma_{1} \end{array}\right]$, $C_{2}=\left[\begin{array}{cc} 0 & 0\\ \sigma_{s} & 0\\ 0 & \sigma_{w} \end{array}\right]$ and $S=\left[\begin{array}{ccc} 0 & 1 & 0\\ 0 & 0 & 1 \end{array}\right]$
b. Let $x_{t}=\left[\begin{array}{c} z_{t}\\ n_{t} \end{array}\right]$ and $u_{t}=\left[\begin{array}{c} n_{t+1}-n_{t}\end{array}\right]$. Then, we can map the problem into
$$\min\mathbb{E}\left[\sum_{t=0}^{\infty}\beta^{t}r\left(x_{t},u_{t}\right)\right]$$where
$$\begin{eqnarray} r\left(x_{t},u_{t}\right) &=& x_{t}'R_{f}x_{t}+u_{t}'Q_{f}u_{t}+2u_{t}'N_{f}x_{t} \\ x_{t+1} &=& Ax_{t}+Bu_{t}+C\epsilon_{t+1} \end{eqnarray}$$by letting $R_{f}=\left[\begin{array}{cccc} 0 & 0 & 0 & 0\\ 0 & 0 & 0 & -\frac{1}{2}\\ 0 & 0 & 0 & \frac{1}{2}\\ 0 & -\frac{1}{2} & \frac{1}{2} & f \end{array}\right]$, $Q_{f}=\left[\begin{array}{c} d\end{array}\right]$, $N_{f}=\left[\begin{array}{cccc} 0 & 0 & 0 & 0\end{array}\right]$, $A=\left[\begin{array}{cc} A_{22} & 0\\ 0 & 1 \end{array}\right]$, $B=\left[\begin{array}{c} 0\\ 1 \end{array}\right]$ and $C=\left[\begin{array}{c} C_{2}\\ 0 \end{array}\right]$.
The optimal policy is given by $u_{t}=-F_{f}x_{t}$. Labor supply is given by $\left[\begin{array}{cccc} 0 & 0 & 0 & 1\end{array}\right]\left(A-BF_{f}\right)x_{t}$.
a. Let $z_{t}=\left[\begin{array}{c} 1\\ b_{t}\\ w_{t} \end{array}\right]$ and $\epsilon_{t}=\left[\begin{array}{c} \epsilon_{1,t}\\ \epsilon_{2,t} \end{array}\right]$. Then, we can obtain
$$\begin{eqnarray} z_{t+1} &=& A_{22}z_{t}+C_{2}\epsilon_{t} \\ \left[\begin{array}{c} b_{t}\\ w_{t} \end{array}\right] &=& Sz_{t} \end{eqnarray}$$by letting $A_{22}=\left[\begin{array}{ccc} 1 & 0 & 0\\ \alpha_{b} & \rho_{1} & 0\\ \alpha_{w} & 0 & \gamma_{1} \end{array}\right]$, $C_{2}=\left[\begin{array}{cc} 0 & 0\\ \sigma_{b} & 0\\ 0 & \sigma_{w} \end{array}\right]$ and $S=\left[\begin{array}{ccc} 0 & 1 & 0\\ 0 & 0 & 1 \end{array}\right]$
b. Let $x_{t}=\left[\begin{array}{c} z_{t}\\ n_{t} \end{array}\right]$ and $u_{t}=\left[\begin{array}{c} n_{t+1}-n_{t}\end{array}\right]$. Then, we can map the problem into
$$\min\mathbb{E}\left[\sum_{t=0}^{\infty}\beta^{t}r\left(x_{t},u_{t}\right)\right]$$where
$$\begin{eqnarray} r\left(x_{t},u_{t}\right) &=& x_{t}'R_{w}x_{t}+u_{t}'Q_{w}u_{t}+2u_{t}'N_{w}x_{t} \\ x_{t+1} &=& Ax_{t}+Bu_{t}+C\epsilon_{t+1} \end{eqnarray}$$by letting $R_{w}=\left[\begin{array}{cccc} 0 & 0 & 0 & 0\\ 0 & 0 & 0 & \frac{1}{2}\\ 0 & 0 & 0 & -\frac{1}{2}\\ 0 & \frac{1}{2} & -\frac{1}{2} & f+4q \end{array}\right]$, $Q_{w}=\left[\begin{array}{c} q\end{array}\right]$, $N_{w}=\left[\begin{array}{cccc} 0 & 0 & 0 & 2q\end{array}\right]$, $A=\left[\begin{array}{cc} A_{22} & 0\\ 0 & 1 \end{array}\right]$, $B=\left[\begin{array}{c} 0\\ 1 \end{array}\right]$, and $C=\left[\begin{array}{c} C_{2}\\ 0 \end{array}\right]$.
The optimal policy is given by $u_{t}=-F_{w}x_{t}$. Labor supply is given by $\left[\begin{array}{cccc} 0 & 0 & 0 & 1\end{array}\right]\left(A-BF_{w}\right)x_{t}$.
Let $x_{t}=\left[\begin{array}{c} z_{t}\\ n_{t} \end{array}\right]$ and $u_{t}=\left[\begin{array}{c} n_{t+1}-n_{t}\end{array}\right]$. Then, we can map the problem into
$$\min\mathbb{E}\left[\sum_{t=0}^{\infty}\beta^{t}r\left(x_{t},u_{t}\right)\right]$$where
$$\begin{eqnarray} r\left(x_{t},u_{t}\right) &=& x_{t}'R_{p}x_{t}+u_{t}'Q_{p}u_{t}+2u_{t}'N_{p}x_{t} \\ x_{t+1} &=& A_{p}x_{t}+B_{p}u_{t}+C_{p}\epsilon_{t+1} \end{eqnarray}$$by letting $R_{p}=\left[\begin{array}{cc} 0 & S'\left[\begin{array}{c} -\frac{1}{2}\\ \frac{1}{2} \end{array}\right]\\ \left[\begin{array}{cc} -\frac{1}{2} & \frac{1}{2}\end{array}\right]S & f+h+4q \end{array}\right]$, $Q_{p}=\left[\begin{array}{c} d+q\end{array}\right]$, $N_{p}=\left[\begin{array}{cccc} 0 & 0 & 0 & 2q\end{array}\right]$, $A_{p}=\left[\begin{array}{cc} A_{22} & 0\\ 0 & 1 \end{array}\right]$, $B_{p}=\left[\begin{array}{c} 0\\ 1 \end{array}\right]$ and $C_{p}=\left[\begin{array}{c} C_{2}\\ 0 \end{array}\right]$.
The optimal policy is given by $u_{t}=-F_{p}x_{t}$. Labor supply is given by $\left[\begin{array}{cccc} 0 & 1\end{array}\right]\left(A_{p}-B_{p}F_{p}\right)x_{t}$.
import numpy as np
import matplotlib.pyplot as plt
import quantecon as qe
np.set_printoptions(precision=4)
# Parameters
β = 0.9
α_s = 5.
ρ_1 = 0.
σ_s = 0.
α_b = 0.
ρ_2 = 0.
σ_b = 0.
f = 1.
h = 1.
q = 0.
d = 0.
# Planner LQ problem
R_tilde = np.array([[0., 0., 0., 0.],
[0., 0., 0., -1/2],
[0., 0., 0., 1/2],
[0., -1/2, 1/2, f + h + 4 * q]])
Q_tilde = np.array([[d + q]])
N_tilde = np.array([[0., 0., 0., 2 * q]])
A_tilde = np.array([[1., 0., 0., 0.],
[α_s, ρ_1, 0., 0.],
[α_b, 0., ρ_2, 0.],
[0., 0., 0., 1.]])
B_tilde = np.array([[0.],
[0.],
[0.],
[1.]])
C_tilde = np.array([[0., 0.],
[σ_s, 0.],
[0., σ_b],
[0., 0.]])
lq_tilde = qe.LQ(Q_tilde, R_tilde, A_tilde, B_tilde, C=C_tilde, N=N_tilde, beta=β)
P_tilde, F_tilde, d_tilde = lq_tilde.stationary_values()
x_0_tilde = np.array([[1., 5., 0., 1.25]])
x_t_tilde = lq_tilde.compute_sequence(x_0_tilde, ts_length=1000)[0]
n_t = x_t_tilde[3]
plt.plot(n_t.ravel(), label='labor')
plt.legend();
P_tilde[[-1]] @ x_t_tilde
array([[0., 0., 0., ..., 0., 0., 0.]])
mask = np.array([1., 1., 0., 1/2])
w_t = -2 * mask * P_tilde[[-1]] @ x_t_tilde
plt.plot(w_t.ravel(), label='wage')
plt.legend();
w_t
array([[2.5, 2.5, 2.5, ..., 2.5, 2.5, 2.5]])
Note:
\boldsymbol{0} & 2q\end{array}\right]$ in the pseudo-code?
LQ
class does minimization so all the signs need to be reversedThe dimensions of your matrices don't match in the Sylvester equation. Assuming $z_t$ has length $n$, the $A_p$ is an $n$ by $n$ matrix and therefore, so is $A_p-B_pF_p$. However, if $\tilde{x}_t'\tilde{R_f}\tilde{x}_t$ is valid, $\tilde{R_f}$ must be $n+1$ by $n+1$. Therefore, $\tilde{R}_{f}+F_{p}'\tilde{Q}_{f}F_{p}+\beta A_{p}'P_{w}\left(A_{p}-B_{p}F_{p}\right)$ must be invalid. I'm guessing what you meant is $z_t'\tilde{R_f}z_t=s_{t}n_{t}-f\tilde{n}_{t}^{2}$. This is achived by setting $\tilde{R}_{f}=\left[\begin{array}{cccc} 0 & 0 & 0 & 0\\ 0 & 0 & 0 & -\frac{1}{2}\\ 0 & 0 & 0 & 0\\ 0 & -\frac{1}{2} & 0 & f \end{array}\right]$ and $\tilde{Q}_{f}=\left[\begin{array}{c} d\end{array}\right]$.
I'm guessing that what you want to do is evaluate the policy for the firm implied by solving the social planner problem. If so, I believe that the equation we want to solve is the following Lyapunov equation:
$$P_{w}=\tilde{R}_{f}+F_{p}'\tilde{Q}_{f}F_{p}+\beta\left(A_{p}-B_{p}F_{p}\right)'P_{w}\left(A_{p}-B_{p}F_{p}\right)$$I'll also solve your Sylvester equation first in case I'm wrong.
R_tilde_f = np.array([[0., 0., 0., 0.],
[0., 0., 0., -1/2],
[0., 0., 0., 0.],
[0., -1/2, 0., f],])
Q_tilde_f = np.array([[d]])
SciPy has a solver for Sylvester equations of the form:
$$AX+XB=Q$$Therefore, we want to find matrices $A_{S}$, $B_{S}$ and $Q_{S}$ such that:
$$A_{S}P_{w}+P_{w}B_{S}=Q_{S}\Longleftrightarrow P_{w}=\tilde{R}_{f}+F_{p}'\tilde{Q}_{f}F_{p}+\beta A_{p}'P_{w}\left(A_{p}-B_{p}F_{p}\right)$$Letting $A_{S}=\beta A_{p}'$, $Q_{S}=-\left(\tilde{R}_{f}+F_{p}'\tilde{Q}_{f}F_{p}\right)\left(A_{p}-B_{p}F_{p}\right)^{-1}$ (assuming the inverse exists) and $B_{S}=-\left(A_{p}-B_{p}F_{p}\right)^{-1}$, we obtain:
$$\beta A_{p}'P_{w}-P_{w}\left(A_{p}-B_{p}F_{p}\right)^{-1}=-\left(\tilde{R}_{f}+F_{p}'\tilde{Q}_{f}F_{p}\right)\left(A_{p}-B_{p}F_{p}\right)^{-1}$$Post-multiplying by $\left(A_{p}-B_{p}F_{p}\right)$, we obtain the required equation.
Alternatively, we can let $A_{S}=-\left(\beta A_{p}'\right)^{-1}$, $B_{S}=\left(A_{p}-B_{p}F_{p}\right)$ and $Q_{S}=-\left(\beta A_{p}'\right)^{-1}\left(\tilde{R}_{f}+F_{p}'\tilde{Q}_{f}F_{p}\right)$.
Note: The inverses don't necessarily exists depending on parameter values.
from scipy.linalg import solve_sylvester
A_sylvester = -np.linalg.inv(β * A_tilde.T)
B_sylvester = A_tilde - B_tilde @ F_tilde
Q_sylvester = A_sylvester @ (R_tilde_f + F_tilde.T @ Q_tilde_f @ F_tilde)
P_sylv = solve_sylvester(A_sylvester, B_sylvester, Q_sylvester)
--------------------------------------------------------------------------- LinAlgError Traceback (most recent call last) <ipython-input-10-905271a33b0f> in <module> 1 from scipy.linalg import solve_sylvester 2 ----> 3 A_sylvester = -np.linalg.inv(β * A_tilde.T) 4 B_sylvester = A_tilde - B_tilde @ F_tilde 5 Q_sylvester = A_sylvester @ (R_tilde_f + F_tilde.T @ Q_tilde_f @ F_tilde) /anaconda3/lib/python3.7/site-packages/numpy/linalg/linalg.py in inv(a) 530 signature = 'D->D' if isComplexType(t) else 'd->d' 531 extobj = get_linalg_error_extobj(_raise_linalgerror_singular) --> 532 ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj) 533 return wrap(ainv.astype(result_t, copy=False)) 534 /anaconda3/lib/python3.7/site-packages/numpy/linalg/linalg.py in _raise_linalgerror_singular(err, flag) 87 88 def _raise_linalgerror_singular(err, flag): ---> 89 raise LinAlgError("Singular matrix") 90 91 def _raise_linalgerror_nonposdef(err, flag): LinAlgError: Singular matrix
P_sylv = np.eye(4)
max_iter = 1000
δ = 1
for i in range(max_iter):
δ = (R_tilde_f + F_tilde.T @ Q_tilde_f @ F_tilde + β * A_tilde.T @ P_sylv @ (A_tilde - B_tilde @ F_tilde)) - P_sylv
P_sylv += δ
if np.max(np.abs(δ)) < 1e-12:
break
if i == max_iter - 1:
print('Max number of iterations reached.')
-2 * P_sylv[[-1]] @ x_t_tilde
array([[25., 25., 25., ..., 25., 25., 25.]])
-2 * P_sylv[[-1]] @ x_t_tilde == w_t
array([[False, False, False, ..., False, False, False]])
This doesn't recover the wage computed manually above.
from quantecon import solve_discrete_lyapunov
A_lyapunov = np.sqrt(β) * (A_tilde - B_tilde @ F_tilde)
B_lyapunov = R_tilde_f + F_tilde.T @ Q_tilde_f @ F_tilde
P_w_2 = solve_discrete_lyapunov(A_lyapunov, B_lyapunov)
-2 * P_w_2[[-1]] @ x_t_tilde
array([[2.5, 2.5, 2.5, ..., 2.5, 2.5, 2.5]])
# Check against the one that I computed by hand
-2 * P_w_2[[-1]] @ x_t_tilde == w_t
array([[ True, True, True, ..., True, True, True]])
Now we look ready to try more complex parametrizations.
TODO: Refactor the code to avoid repetition
# Parameters
β = 0.9
α_s = 5.
ρ_1 = 0.
σ_s = 0.
α_b = 0.
ρ_2 = 0.
σ_b = 0.
f = 1.
h = 1.
q = 0.
d = 1.
# Planner LQ problem
R_tilde = np.array([[0., 0., 0., 0.],
[0., 0., 0., -1/2],
[0., 0., 0., 1/2],
[0., -1/2, 1/2, f + h + 4 * q]])
Q_tilde = np.array([[d + q]])
N_tilde = np.array([[0., 0., 0., 2 * q]])
A_tilde = np.array([[1., 0., 0., 0.],
[α_s, ρ_1, 0., 0.],
[α_b, 0., ρ_2, 0.],
[0., 0., 0., 1.]])
B_tilde = np.array([[0.],
[0.],
[0.],
[1.]])
C_tilde = np.array([[0., 0.],
[σ_s, 0.],
[0., σ_b],
[0., 0.]])
lq_tilde = qe.LQ(Q_tilde, R_tilde, A_tilde, B_tilde, C=C_tilde, N=N_tilde, beta=β)
P_tilde, F_tilde, d_tilde = lq_tilde.stationary_values()
Check the following:
$$\tilde{P}=\tilde{R}+\tilde{F}'\tilde{Q}\tilde{F}+\beta\left(\tilde{A}-\tilde{B}\tilde{F}\right)'\tilde{P}\left(\tilde{A}-\tilde{B}\tilde{F}\right)$$flow_utility = R_tilde + F_tilde.T @ Q_tilde @ F_tilde - 2 * F_tilde.T @ N_tilde
continuation_val = (A_tilde - B_tilde @ F_tilde).T @ P_tilde @ (A_tilde - B_tilde @ F_tilde)
# Should be very close to 0 if above holds
np.max(np.abs(flow_utility + β * continuation_val - P_tilde))
2.220446049250313e-16
# Should be very close to 0
np.max(np.abs(solve_discrete_lyapunov(A=np.sqrt(β) * (A_tilde - B_tilde @ F_tilde).T,
B=R_tilde + F_tilde.T @ Q_tilde @ F_tilde - 2 * F_tilde.T @ N_tilde)
- P_tilde))
5.684341886080802e-14
# Break down the problem for the firm
R_tilde_f = np.array([[0., 0., 0., 0.],
[0., 0., 0., -1/2],
[0., 0., 0., 0.],
[0., -1/2, 0., f],])
Q_tilde_f = np.array([[d]])
A_lyapunov = np.sqrt(β) * (A_tilde - B_tilde @ F_tilde).T
B_lyapunov = R_tilde_f + F_tilde.T @ Q_tilde_f @ F_tilde
P_w = solve_discrete_lyapunov(A=A_lyapunov, B=B_lyapunov)
# Same approach for the worker
R_tilde_w = np.array([[0., 0., 0., 0.],
[0., 0., 0., 0.],
[0., 0., 0., 1/2],
[0., 0., 1/2, h + 4 * q]])
Q_tilde_w = np.array([[q]])
N_tilde_w = np.array([[0., 0., 0., 2 * q]])
B_lyapunov2 = R_tilde_w + F_tilde.T @ Q_tilde_w @ F_tilde - 2 * F_tilde.T @ N_tilde_w
P_w2 = solve_discrete_lyapunov(A=A_lyapunov, B=B_lyapunov2)
R_tilde_w + R_tilde_f == R_tilde
array([[ True, True, True, True], [ True, True, True, True], [ True, True, True, True], [ True, True, True, True]])
Q_tilde_w + Q_tilde_f == Q_tilde
array([[ True]])
N_tilde_w == N_tilde
array([[ True, True, True, True]])
# Should be very close to 0
P_w3 = solve_discrete_lyapunov(A=A_lyapunov, B=B_lyapunov+B_lyapunov2)
np.max(np.abs(P_w3 - P_tilde))
5.684341886080802e-14
Implied wages
- 2 * P_w[[-1]] @ x_0_tilde.T
array([[3.3864]])
2 * P_w2[[-1]] @ x_0_tilde.T
array([[3.3864]])
2 * P_w3[[-1]] @ x_0_tilde.T
array([[0.]])
Do value functions add up?
x_0_tilde @ P_w @ x_0_tilde.T + x_0_tilde @ P_w2 @ x_0_tilde.T
array([[-31.25]])
x_0_tilde @ P_w3 @ x_0_tilde.T
array([[-31.25]])
# Wage function with d=1 (computed by hand)
(np.array([0., 1., 0., -2]) - 2 * F_tilde)
array([[ 1.7729, 1. , 0. , -3.4183]])
(np.array([0., 1., 0., -2]) - 2 * F_tilde) @ x_0_tilde.T
array([[2.5]])
# Firm LQ problem
S_w = -2 * P_w[[-1]]
S_s = np.array([[0., 1., 0., 0.]])
cross_term = -(S_s - S_w) / 2
n = A_tilde.shape[0]
R_hat = np.block([[np.zeros((n, n)), cross_term.T],
[cross_term, f]])
Q_hat = np.array([[d]])
N_hat = np.zeros((1, n + 1))
A_hat = np.block([[A_tilde - B_tilde @ F_tilde, np.zeros((n, 1))],
[np.zeros((1, n)), np.ones((1, 1))]])
B_hat = np.array([[0.],
[0.],
[0.],
[0.],
[1.]])
C_hat = np.vstack((C_tilde, np.zeros((1, C_tilde.shape[1]))))
lq_hat = qe.LQ(Q_hat, R_hat, A_hat, B_hat, C=C_hat, N=N_hat, beta=β)
P_hat, F_hat, d_hat = lq_hat.stationary_values()
F_tilde
array([[-0.8864, 0. , 0. , 0.7092]])
F_hat
array([[-0.2292, 0. , 0. , -0.1964, 0.5884]])
x_0_hat = np.array([[1., 5., 0., 1.25, 1.25]]).T
x_t_hat = lq_hat.compute_sequence(x_0_hat, ts_length=1000)[0]
plt.plot(x_t_hat[3], label='planner implied')
plt.plot(x_t_hat[4], label='firm implied');
plt.legend();
x_t_hat[3]
array([1.25, 1.25, 1.25, ..., 1.25, 1.25, 1.25])
np.var(x_t_hat[3])
0.0
x_t_hat[4]
array([1.25 , 0.9892, 0.8819, ..., 0.8068, 0.8068, 0.8068])
np.var(x_t_hat[4])
0.0002357102288177184
# Sup norm
np.max(np.abs(x_t_hat[3] - x_t_hat[4]))
0.44322084543531903
# Worker LQ problem
S_b = np.array([[0., 0., 1., 0.]])
cross_term = (S_b - S_w) / 2
n = A_tilde.shape[0]
R_hat = np.block([[np.zeros((n, n)), cross_term.T],
[cross_term, h + 4 * q]])
Q_hat = np.array([[q]])
N_hat = np.array([[0., 0., 0., 0., 2 * q]])
A_hat = np.block([[A_tilde - B_tilde @ F_tilde, np.zeros((n, 1))],
[np.zeros((1, n)), np.ones((1, 1))]])
B_hat = np.array([[0.],
[0.],
[0.],
[0.],
[1.]])
C_hat = np.vstack((C_tilde, np.zeros((1, C_tilde.shape[1]))))
lq_hat = qe.LQ(Q_hat, R_hat, A_hat, B_hat, C=C_hat, N=N_hat, beta=β)
P_hat, F_hat, d_hat = lq_hat.stationary_values()
x_0_hat = np.array([[1., 5., 0., 1.25, 1.25]]).T
x_t_hat = lq_hat.compute_sequence(x_0_hat, ts_length=1000)[0]
plt.plot(x_t_hat[3], label='planner implied')
plt.plot(x_t_hat[4], label='worked implied')
plt.legend();
# Sup norm
np.max(np.abs(x_t_hat[3] - x_t_hat[4]))
0.44322084543531926