hw4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 1\n",
    "\n",
    "**Note:** \n",
    "- $w_t$ is used as a state variable in the problem but as a shock variable in the QuantEcon documentation. This might potentially confuse some students.\n",
    "- In exercise 3, you use $A_p$, $B_p$ and $C_p$ as notation -- did you mean to use $A_f$, $B_f$ and $C_f$ here?\n",
    "\n",
    "**a.** Let $z_{t}=\\left[\\begin{array}{c}\n",
    "1\\\\\n",
    "s_{t}\\\\\n",
    "w_{t}\n",
    "\\end{array}\\right]$ and $\\epsilon_{t}=\\left[\\begin{array}{c}\n",
    "\\epsilon_{1,t}\\\\\n",
    "\\epsilon_{2,t}\n",
    "\\end{array}\\right]$. Then, we can obtain \n",
    "\n",
    "$$\\begin{eqnarray}\n",
    "z_{t+1} &=& A_{22}z_{t}+C_{2}\\epsilon_{t} \\\\\n",
    "\\left[\\begin{array}{c}\n",
    "s_{t}\\\\\n",
    "w_{t}\n",
    "\\end{array}\\right] &=& Sz_{t}\n",
    "\\end{eqnarray}$$\n",
    "\n",
    "by letting $A_{22}=\\left[\\begin{array}{ccc}\n",
    "1 & 0 & 0\\\\\n",
    "\\alpha_{s} & \\rho_{1} & 0\\\\\n",
    "\\alpha_{w} & 0 & \\gamma_{1}\n",
    "\\end{array}\\right]$, $C_{2}=\\left[\\begin{array}{cc}\n",
    "0 & 0\\\\\n",
    "\\sigma_{s} & 0\\\\\n",
    "0 & \\sigma_{w}\n",
    "\\end{array}\\right]$ and $S=\\left[\\begin{array}{ccc}\n",
    "0 & 1 & 0\\\\\n",
    "0 & 0 & 1\n",
    "\\end{array}\\right]$\n",
    "\n",
    "**b.** Let $x_{t}=\\left[\\begin{array}{c}\n",
    "z_{t}\\\\\n",
    "n_{t}\n",
    "\\end{array}\\right]$ and $u_{t}=\\left[\\begin{array}{c}\n",
    "n_{t+1}-n_{t}\\end{array}\\right]$. Then, we can map the problem into \n",
    "\n",
    "$$\\min\\mathbb{E}\\left[\\sum_{t=0}^{\\infty}\\beta^{t}r\\left(x_{t},u_{t}\\right)\\right]$$ \n",
    "\n",
    "where \n",
    "\n",
    "$$\\begin{eqnarray}\n",
    "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} \\\\ \n",
    "x_{t+1} &=& Ax_{t}+Bu_{t}+C\\epsilon_{t+1}\n",
    "\\end{eqnarray}$$\n",
    "\n",
    "by letting $R_{f}=\\left[\\begin{array}{cccc}\n",
    "0 & 0 & 0 & 0\\\\\n",
    "0 & 0 & 0 & -\\frac{1}{2}\\\\\n",
    "0 & 0 & 0 & \\frac{1}{2}\\\\\n",
    "0 & -\\frac{1}{2} & \\frac{1}{2} & f\n",
    "\\end{array}\\right]$, $Q_{f}=\\left[\\begin{array}{c}\n",
    "d\\end{array}\\right]$, $N_{f}=\\left[\\begin{array}{cccc}\n",
    "0 & 0 & 0 & 0\\end{array}\\right]$, $A=\\left[\\begin{array}{cc}\n",
    "A_{22} & 0\\\\\n",
    "0 & 1\n",
    "\\end{array}\\right]$, $B=\\left[\\begin{array}{c}\n",
    "0\\\\\n",
    "1\n",
    "\\end{array}\\right]$ and $C=\\left[\\begin{array}{c}\n",
    "C_{2}\\\\\n",
    "0\n",
    "\\end{array}\\right]$.\n",
    "\n",
    "The optimal policy is given by $u_{t}=-F_{f}x_{t}$.\n",
    "Labor supply is given by $\\left[\\begin{array}{cccc}\n",
    "0 & 0 & 0 & 1\\end{array}\\right]\\left(A-BF_{f}\\right)x_{t}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2\n",
    "\n",
    "**a.** Let $z_{t}=\\left[\\begin{array}{c}\n",
    "1\\\\\n",
    "b_{t}\\\\\n",
    "w_{t}\n",
    "\\end{array}\\right]$ and $\\epsilon_{t}=\\left[\\begin{array}{c}\n",
    "\\epsilon_{1,t}\\\\\n",
    "\\epsilon_{2,t}\n",
    "\\end{array}\\right]$. Then, we can obtain \n",
    "\n",
    "$$\\begin{eqnarray}\n",
    "z_{t+1} &=& A_{22}z_{t}+C_{2}\\epsilon_{t} \\\\\n",
    "\\left[\\begin{array}{c}\n",
    "b_{t}\\\\\n",
    "w_{t}\n",
    "\\end{array}\\right] &=& Sz_{t}\n",
    "\\end{eqnarray}$$\n",
    "\n",
    "by letting $A_{22}=\\left[\\begin{array}{ccc}\n",
    "1 & 0 & 0\\\\\n",
    "\\alpha_{b} & \\rho_{1} & 0\\\\\n",
    "\\alpha_{w} & 0 & \\gamma_{1}\n",
    "\\end{array}\\right]$, $C_{2}=\\left[\\begin{array}{cc}\n",
    "0 & 0\\\\\n",
    "\\sigma_{b} & 0\\\\\n",
    "0 & \\sigma_{w}\n",
    "\\end{array}\\right]$ and $S=\\left[\\begin{array}{ccc}\n",
    "0 & 1 & 0\\\\\n",
    "0 & 0 & 1\n",
    "\\end{array}\\right]$\n",
    "\n",
    "**b.** Let $x_{t}=\\left[\\begin{array}{c}\n",
    "z_{t}\\\\\n",
    "n_{t}\n",
    "\\end{array}\\right]$ and $u_{t}=\\left[\\begin{array}{c}\n",
    "n_{t+1}-n_{t}\\end{array}\\right]$. Then, we can map the problem into \n",
    "\n",
    "$$\\min\\mathbb{E}\\left[\\sum_{t=0}^{\\infty}\\beta^{t}r\\left(x_{t},u_{t}\\right)\\right]$$ \n",
    "\n",
    "where \n",
    "\n",
    "$$\\begin{eqnarray}\n",
    "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} \\\\ \n",
    "x_{t+1} &=& Ax_{t}+Bu_{t}+C\\epsilon_{t+1}\n",
    "\\end{eqnarray}$$\n",
    "\n",
    "by letting $R_{w}=\\left[\\begin{array}{cccc}\n",
    "0 & 0 & 0 & 0\\\\\n",
    "0 & 0 & 0 & \\frac{1}{2}\\\\\n",
    "0 & 0 & 0 & -\\frac{1}{2}\\\\\n",
    "0 & \\frac{1}{2} & -\\frac{1}{2} & f+4q\n",
    "\\end{array}\\right]$, $Q_{w}=\\left[\\begin{array}{c}\n",
    "q\\end{array}\\right]$, $N_{w}=\\left[\\begin{array}{cccc}\n",
    "0 & 0 & 0 & 2q\\end{array}\\right]$, $A=\\left[\\begin{array}{cc}\n",
    "A_{22} & 0\\\\\n",
    "0 & 1\n",
    "\\end{array}\\right]$, $B=\\left[\\begin{array}{c}\n",
    "0\\\\\n",
    "1\n",
    "\\end{array}\\right]$, and $C=\\left[\\begin{array}{c}\n",
    "C_{2}\\\\\n",
    "0\n",
    "\\end{array}\\right]$.\n",
    "\n",
    "The optimal policy is given by $u_{t}=-F_{w}x_{t}$.\n",
    "Labor supply is given by $\\left[\\begin{array}{cccc}\n",
    "0 & 0 & 0 & 1\\end{array}\\right]\\left(A-BF_{w}\\right)x_{t}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3\n",
    "\n",
    "Let $x_{t}=\\left[\\begin{array}{c}\n",
    "z_{t}\\\\\n",
    "n_{t}\n",
    "\\end{array}\\right]$ and $u_{t}=\\left[\\begin{array}{c}\n",
    "n_{t+1}-n_{t}\\end{array}\\right]$. Then, we can map the problem into \n",
    "\n",
    "$$\\min\\mathbb{E}\\left[\\sum_{t=0}^{\\infty}\\beta^{t}r\\left(x_{t},u_{t}\\right)\\right]$$ \n",
    "\n",
    "where \n",
    "\n",
    "$$\\begin{eqnarray}\n",
    "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} \\\\ \n",
    "x_{t+1} &=& A_{p}x_{t}+B_{p}u_{t}+C_{p}\\epsilon_{t+1}\n",
    "\\end{eqnarray}$$\n",
    "\n",
    "by letting $R_{p}=\\left[\\begin{array}{cc}\n",
    "0 & S'\\left[\\begin{array}{c}\n",
    "-\\frac{1}{2}\\\\\n",
    "\\frac{1}{2}\n",
    "\\end{array}\\right]\\\\\n",
    "\\left[\\begin{array}{cc}\n",
    "-\\frac{1}{2} & \\frac{1}{2}\\end{array}\\right]S & f+h+4q\n",
    "\\end{array}\\right]$, $Q_{p}=\\left[\\begin{array}{c}\n",
    "d+q\\end{array}\\right]$, $N_{p}=\\left[\\begin{array}{cccc}\n",
    "0 & 0 & 0 & 2q\\end{array}\\right]$, $A_{p}=\\left[\\begin{array}{cc}\n",
    "A_{22} & 0\\\\\n",
    "0 & 1\n",
    "\\end{array}\\right]$, $B_{p}=\\left[\\begin{array}{c}\n",
    "0\\\\\n",
    "1\n",
    "\\end{array}\\right]$ and $C_{p}=\\left[\\begin{array}{c}\n",
    "C_{2}\\\\\n",
    "0\n",
    "\\end{array}\\right]$.\n",
    "\n",
    "The optimal policy is given by $u_{t}=-F_{p}x_{t}$.\n",
    "Labor supply is given by $\\left[\\begin{array}{cccc}\n",
    "0 & 1\\end{array}\\right]\\left(A_{p}-B_{p}F_{p}\\right)x_{t}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import quantecon as qe\n",
    "\n",
    "np.set_printoptions(precision=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "β = 0.9\n",
    "\n",
    "α_s = 5.\n",
    "ρ_1 = 0.\n",
    "σ_s = 0.\n",
    "\n",
    "α_b = 0.\n",
    "ρ_2 = 0.\n",
    "σ_b = 0.\n",
    "\n",
    "f = 1.\n",
    "h = 1.\n",
    "q = 0.\n",
    "d = 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Planner LQ problem\n",
    "R_tilde = np.array([[0., 0., 0., 0.],\n",
    "                    [0., 0., 0., -1/2],\n",
    "                    [0., 0., 0., 1/2],\n",
    "                    [0., -1/2, 1/2, f + h + 4 * q]])\n",
    "\n",
    "Q_tilde = np.array([[d + q]])\n",
    "\n",
    "N_tilde = np.array([[0., 0., 0., 2 * q]])\n",
    "\n",
    "A_tilde = np.array([[1., 0., 0., 0.],\n",
    "                    [α_s, ρ_1, 0., 0.],\n",
    "                    [α_b, 0., ρ_2, 0.],\n",
    "                    [0., 0., 0., 1.]])\n",
    "\n",
    "B_tilde = np.array([[0.],\n",
    "                    [0.],\n",
    "                    [0.],\n",
    "                    [1.]])\n",
    "\n",
    "C_tilde = np.array([[0., 0.],\n",
    "                    [σ_s, 0.],\n",
    "                    [0., σ_b],\n",
    "                    [0., 0.]])\n",
    "\n",
    "lq_tilde = qe.LQ(Q_tilde, R_tilde, A_tilde, B_tilde, C=C_tilde, N=N_tilde, beta=β)\n",
    "P_tilde, F_tilde, d_tilde = lq_tilde.stationary_values()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_0_tilde = np.array([[1., 5., 0., 1.25]])\n",
    "\n",
    "x_t_tilde = lq_tilde.compute_sequence(x_0_tilde, ts_length=1000)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAEkpJREFUeJzt3X+MndV95/H3J3iINxAasCerBOMdZxtqKAwFXaWkWAvUNKW/YKsSJYgfTiFCiZptqCo1P1bCSvhrpapNVypQi19RipxVKGlSdikbuUgIObAZl4gaHBO3ODCQxoPNkqVZEpx+94957Mw69sx45prL3PN+SVe+zznnPvd75ow/89znPncmVYUkqR1vGnQBkqTXl8EvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JasyyQRdwOCtXrqyxsbFBlyFJS8a2bdterKrR+Yx9Qwb/2NgYExMTgy5DkpaMJN+Z71hP9UhSYwx+SWqMwS9JjXlDnuOXpPl67bXXmJyc5NVXXx10Ka+L5cuXs2rVKkZGRha8D4Nf0pI2OTnJW9/6VsbGxkgy6HKOqapi7969TE5OsmbNmgXvx1M9kpa0V199lRUrVgx96AMkYcWKFYt+dWPwS1ryWgj9A/oxV4Nfkhpj8EvSIp144omz9u/evZuzzjrrdapmbga/JL3B/fjHP+7r/gx+SeqTV155hfXr13Peeedx9tln85WvfOVg3/79+9mwYQPj4+NcccUV/OAHPwBgy5YtnHvuuZx99tlcd911/PCHPwSmf3XNZz/7WdatW8eXvvSlvtbp5ZyShsZn/uZJnnrh+33d55nvPImNv/Xz8xq7fPlyvvzlL3PSSSfx4osvcv7553PZZZcBsHPnTu644w4uuOACrrvuOm655RY+9rGP8aEPfYgtW7Zw+umnc+2113Lrrbdy4403HtzfI4880tf5gEf8ktQ3VcWnP/1pxsfHueSSS3j++ef53ve+B8Bpp53GBRdcAMDVV1/NI488ws6dO1mzZg2nn346ABs2bODhhx8+uL8PfOADx6ROj/glDY35HpkfK/fccw9TU1Ns27aNkZERxsbGDl5zf+hlmEmoqln3d8IJJxyTOj3il6Q+efnll3n729/OyMgIDz30EN/5zk9+U/Kzzz7L17/+dQA2b97MunXrWLt2Lbt372bXrl0AfOELX+DCCy885nUa/JLUJ1dddRUTExP0ej3uuece1q5de7DvjDPO4POf/zzj4+Ps27ePj370oyxfvpy77rqL97///Zx99tm86U1v4iMf+cgxr9NTPZK0SK+88goAK1euPHhUf6innnrqsO3r16/n8ccf/6n23bt3962+Q3nEL0mNMfglqTEGv6Qlb66rY4ZJP+Zq8Eta0pYvX87evXubCP8Dv49/+fLli9qPb+5KWtJWrVrF5OQkU1NTgy7ldXHgL3AthsEvaUkbGRlZ1F+japGneiSpMfMK/iR3JtmTZPsR+i9P8kSSbyaZSLJuRt+GJN/ubhv6VbgkaWHme8R/N3DpLP1bgHOq6heA64DbAZKcAmwEfhF4D7AxyckLrlaStGjzCv6qehjYN0v/K/WTt9RPAA7c/1Xga1W1r6peAr7G7D9AJEnHWN/O8Sf57STfAv4700f9AKcCz80YNtm1SZIGpG/BX1Vfrqq1wH8Ebu6aD/fn4A97sW2SG7r3ByZauSxLkgah71f1dKeF/n2SlUwf4Z82o3sV8MIRHrepqnpV1RsdHe13WZKkTl+CP8nPpvsrA0nOA44H9gIPAu9LcnL3pu77ujZJ0oDM6wNcSTYDFwErk0wyfaXOCEBV3Qb8DnBtkteA/wt8oHuzd1+Sm4FvdLv6bFUd8U1iSdKxlzfi77fo9Xo1MTEx6DIkaclIsq2qevMZ6yd3JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWrMnMGf5M4ke5JsP0L/VUme6G5bk5wzo+8PkjyZZHuSzUmW97N4SdLRm88R/93ApbP0PwNcWFXjwM3AJoAkpwK/D/Sq6izgOOCDi6pWkrRoy+YaUFUPJxmbpX/rjM1HgVWH7P/fJHkNeAvwwsLKlCT1S7/P8V8PPABQVc8Dfww8C3wXeLmq/mefn0+SdJT6FvxJLmY6+D/RbZ8MXA6sAd4JnJDk6lkef0OSiSQTU1NT/SpLknSIvgR/knHgduDyqtrbNV8CPFNVU1X1GnAf8EtH2kdVbaqqXlX1RkdH+1GWJOkwFh38SVYzHerXVNXTM7qeBc5P8pYkAdYDOxb7fJKkxZnzzd0km4GLgJVJJoGNwAhAVd0G3ASsAG6Zznf2d0fujyW5F/h7YD/wON0VP5KkwUlVDbqGn9Lr9WpiYmLQZUjSkpFkW1X15jPWT+5KUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1Jg5gz/JnUn2JNl+hP6rkjzR3bYmOWdG39uS3JvkW0l2JHlvP4uXJB29+Rzx3w1cOkv/M8CFVTUO3AxsmtH3Z8DfVtVa4BxgxwLrlCT1ybK5BlTVw0nGZunfOmPzUWAVQJKTgP8AfKgb9yPgRwsvVZLUD/0+x3898EB3/13AFHBXkseT3J7khCM9MMkNSSaSTExNTfW5LEnSAX0L/iQXMx38n+ialgHnAbdW1bnAvwCfPNLjq2pTVfWqqjc6OtqvsiRJh+hL8CcZB24HLq+qvV3zJDBZVY912/cy/YNAkjRAiw7+JKuB+4BrqurpA+1V9c/Ac0l+rmtaDzy12OeTJC3OnG/uJtkMXASsTDIJbARGAKrqNuAmYAVwSxKA/VXV6x7+n4B7khwP/BPwu/2egCTp6Mznqp4r5+j/MPDhI/R9E+gdrk+SNBh+cleSGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhoz5yd3l5LP/M2TPPXC9wddhiQtyJnvPImNv/Xzx/x5POKXpMYM1RH/6/GTUpKWOo/4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDVmzuBPcmeSPUm2H6H/qiRPdLetSc45pP+4JI8nub9fRUuSFm4+R/x3A5fO0v8McGFVjQM3A5sO6f84sGNB1UmS+m7O4K+qh4F9s/RvraqXus1HgVUH+pKsAn4DuH2RdUqS+qTf5/ivBx6Ysf054I+Af+3z80iSFqhvwZ/kYqaD/xPd9m8Ce6pq2zwff0OSiSQTU1NT/SpLknSIvgR/knGmT+dcXlV7u+YLgMuS7Aa+CPxykr880j6qalNV9aqqNzo62o+yJEmHsejgT7IauA+4pqqePtBeVZ+qqlVVNQZ8EPi7qrp6sc8nSVqcOf8CV5LNwEXAyiSTwEZgBKCqbgNuAlYAtyQB2F9VvWNVsCRpcVJVg67hp/R6vZqYmBh0GZK0ZCTZNt+Dbj+5K0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1Jj5gz+JHcm2ZNk+xH6r0ryRHfbmuScrv20JA8l2ZHkySQf73fxkqSjN58j/ruBS2fpfwa4sKrGgZuBTV37fuAPq+oM4Hzg95KcuYhaJUl9MGfwV9XDwL5Z+rdW1Uvd5qPAqq79u1X19939/wPsAE5ddMWSpEXp9zn+64EHDm1MMgacCzzW5+eTJB2lZf3aUZKLmQ7+dYe0nwj8FXBjVX1/lsffANwAsHr16n6VJUk6RF+O+JOMA7cDl1fV3hntI0yH/j1Vdd9s+6iqTVXVq6re6OhoP8qSJB3GooM/yWrgPuCaqnp6RnuAO4AdVfUni30eSVJ/zHmqJ8lm4CJgZZJJYCMwAlBVtwE3ASuAW6aznv1V1QMuAK4B/iHJN7vdfbqq/ke/JyFJmr85g7+qrpyj/8PAhw/T/giQhZcmSToW/OSuJDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqzJzBn+TOJHuSbD9C/1VJnuhuW5OcM6Pv0iQ7k+xK8sl+Fi5JWpj5HPHfDVw6S/8zwIVVNQ7cDGwCSHIc8OfArwFnAlcmOXNR1UqSFm3O4K+qh4F9s/RvraqXus1HgVXd/fcAu6rqn6rqR8AXgcsXWa8kaZH6fY7/euCB7v6pwHMz+ia7NknSAC3r146SXMx08K870HSYYTXL428AbgBYvXp1v8qSJB2iL0f8ScaB24HLq2pv1zwJnDZj2CrghSPto6o2VVWvqnqjo6P9KEuSdBiLDv4kq4H7gGuq6ukZXd8A3p1kTZLjgQ8CX13s80mSFmfOUz1JNgMXASuTTAIbgRGAqroNuAlYAdySBGB/d+S+P8nHgAeB44A7q+rJYzILSdK8peqIp90Hptfr1cTExKDLkKQlI8m2qurNZ6yf3JWkxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4JekxqSqBl3DT0kyBXxngQ9fCbzYx3KWAuc8/FqbLzjno/Xvqmp0PgPfkMG/GEkmqqo36DpeT855+LU2X3DOx5KneiSpMQa/JDVmGIN/06ALGADnPPxamy8452Nm6M7xS5JmN4xH/JKkWQxN8Ce5NMnOJLuSfHLQ9fRLktOSPJRkR5Ink3y8az8lydeSfLv79+SuPUn+a/d1eCLJeYOdwcIlOS7J40nu77bXJHmsm/N/S3J81/7mbntX1z82yLoXKsnbktyb5Fvder932Nc5yR9039fbk2xOsnzY1jnJnUn2JNk+o+2o1zXJhm78t5NsWExNQxH8SY4D/hz4NeBM4MokZw62qr7ZD/xhVZ0BnA/8Xje3TwJbqurdwJZuG6a/Bu/ubjcAt77+JffNx4EdM7b/C/Cn3ZxfAq7v2q8HXqqqnwX+tBu3FP0Z8LdVtRY4h+m5D+06JzkV+H2gV1VnAccBH2T41vlu4NJD2o5qXZOcAmwEfhF4D7DxwA+LBamqJX8D3gs8OGP7U8CnBl3XMZrrV4BfAXYC7+ja3gHs7O7/BXDljPEHxy2lG7Cq+w/xy8D9QJj+YMuyQ9cceBB4b3d/WTcug57DUc73JOCZQ+se5nUGTgWeA07p1u1+4FeHcZ2BMWD7QtcVuBL4ixnt/9+4o70NxRE/P/kGOmCyaxsq3Uvbc4HHgH9bVd8F6P59ezdsWL4WnwP+CPjXbnsF8L+ran+3PXNeB+fc9b/cjV9K3gVMAXd1p7duT3ICQ7zOVfU88MfAs8B3mV63bQz3Oh9wtOva1/UeluDPYdqG6nKlJCcCfwXcWFXfn23oYdqW1NciyW8Ce6pq28zmwwytefQtFcuA84Bbq+pc4F/4ycv/w1nyc+5OVVwOrAHeCZzA9KmOQw3TOs/lSHPs69yHJfgngdNmbK8CXhhQLX2XZITp0L+nqu7rmr+X5B1d/zuAPV37MHwtLgAuS7Ib+CLTp3s+B7wtybJuzMx5HZxz1/8zwL7Xs+A+mAQmq+qxbvtepn8QDPM6XwI8U1VTVfUacB/wSwz3Oh9wtOva1/UeluD/BvDu7mqA45l+g+irA66pL5IEuAPYUVV/MqPrq8CBd/Y3MH3u/0D7td3VAecDLx94SblUVNWnqmpVVY0xvZZ/V1VXAQ8BV3TDDp3zga/FFd34JXUkWFX/DDyX5Oe6pvXAUwzxOjN9iuf8JG/pvs8PzHlo13mGo13XB4H3JTm5e6X0vq5tYQb9pkcf3zz5deBp4B+B/zzoevo4r3VMv6R7Avhmd/t1ps9tbgG+3f17Sjc+TF/h9I/APzB9xcTA57GI+V8E3N/dfxfwv4BdwJeAN3fty7vtXV3/uwZd9wLn+gvARLfWfw2cPOzrDHwG+BawHfgC8OZhW2dgM9PvYbzG9JH79QtZV+C6bu67gN9dTE1+cleSGjMsp3okSfNk8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1Jj/ByG9W70aqoUIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_t = x_t_tilde[3]\n",
    "plt.plot(n_t.ravel(), label='labor')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0., 0., 0., ..., 0., 0., 0.]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P_tilde[[-1]] @ x_t_tilde"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAEvNJREFUeJzt3X+QnVWd5/H3Z0kkTEDkR2Q0IXaoRRnjQoINGkMoBJUsxaw7VVgFuzKwWBWKdQZIsbWjbDGUu+WWWsi4OjKAA1K1xQ66wgBDoUgxcYHaHZZOzABJgwLDQC9ZQVSClHES/O4f/STT23TI7R9J033er6pbfZ/zfO/T5/TpfO5zz31uOlWFJKkd/2S6OyBJ2rcMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1Jj5kx3B8Zy+OGHV19f33R3Q5JmjPXr1/+0qhb0UvumDP6+vj4GBgamuxuSNGMk+ftea13qkaTGGPyS1BiDX5Ia86Zc45ekXmzfvp2hoSG2bds23V3ZZ+bNm8eiRYuYO3fuhI9h8EuasYaGhjjooIPo6+sjyXR3Z6+rKl566SWGhoZYsmTJhI/jUo+kGWvbtm0cdthhTYQ+QBIOO+ywSb/CMfglzWithP5OUzFeg1+SGmPwS1JjDH5JaozBL0kT9KUvfYmvfvWrAKxdu5ZTTz0VgPvuu49PfvKTXHTRRfT397N06VKuvPLKXY+7++67OeaYYzjppJO4+OKLOfPMMwF49dVXueCCCzjhhBNYvnw5d9xxx17pt5dzSpoVPvdXm9j8/NYpPeZ73/lWrvzdpbvdf/LJJ/PlL3+Ziy++mIGBAX7961+zfft2HnzwQVatWsUnPvEJDj30UF577TVOO+00HnnkEd797ndz4YUXcv/997NkyRLOOeecXcf7/Oc/z6mnnsqNN97IL37xC0488UQ+8pGPMH/+/Ckdl2f8kjRB73//+1m/fj2vvPIK+++/PytWrGBgYIAHHniAVatW8e1vf5vjjz+e5cuXs2nTJjZv3szjjz/OUUcdtes6/JHB//3vf58vfOELLFu2jFNOOYVt27bx7LPPTnm/PeOXNCu80Zn53jJ37lz6+vr45je/yYc+9CGOPfZY1q1bx1NPPcUBBxzAVVddxcMPP8whhxzC+eefz7Zt26iq3R6vqrj11lt5z3ves1f77Rm/JE3CySefzFVXXcXJJ5/MqlWruPbaa1m2bBlbt25l/vz5HHzwwfzkJz/hu9/9LgDHHHMMTz/9NM888wwA3/rWt3Yd6/TTT+drX/varieHH/7wh3ulzwa/JE3CqlWr2LJlCytWrOCII45g3rx5rFq1iuOOO47ly5ezdOlSLrjgAlauXAnAAQccwDXXXMPq1as56aSTOOKIIzj44IMBuOKKK9i+fTvHHnss73vf+7jiiiv2Sp9d6pGkSTjttNPYvn37ru0f/ehHu+7fdNNNYz7mwx/+MI8//jhVxac//Wn6+/uB4SeF6667bq/2Fzzjl6R97hvf+AbLli1j6dKlvPzyy1x44YX79Pt7xi9J+9jatWtZu3bttH1/z/glzWhvdJXMbDQV4zX4Jc1Y8+bN46WXXmom/Hf+f/zz5s2b1HFc6pE0Yy1atIihoSFefPHF6e7KPrPzL3BNhsEvacaaO3fupP4SVatc6pGkxuwx+JMcmWRdksEkm5Jcspu6U5Js7Gr+x4j21UmeSPJkks9MZeclSePXy1LPDuCyqtqQ5CBgfZJ7q2rzzoIkbwOuAVZX1bNJ3t617wd8HfgoMAQ8nOTOkY+VJO1bezzjr6otVbWhu/8KMAgsHFX2r4DbqurZru6Frv1E4Mmqerqq/gG4Bfj4VHVekjR+41rjT9IHLAceGrXr3cAhSX6QZH2S3+/aFwLPjagb4vVPGpKkfajnq3qSHAjcClxaVaP/2sEc4P3AacABwP9K8jfAWH8OfswLbpOsAdYALF68uNduSZLGqacz/iRzGQ79m6vqtjFKhoDvVdWrVfVT4H7guK79yBF1i4Dnx/oeVXV9VfVXVf+CBQvGMwZJ0jj0clVPgBuAwaq6ejdldwCrksxJ8lvABxh+L+Bh4OgkS5K8BTgbuHNqui5JmohelnpWAucCjybZ2LVdDiwGqKprq2owyfeAR4DfAH9eVY8BJPkD4B5gP+DGqto0xWOQJI1D3oz/x0V/f38NDAxMdzckacZIsr6q+nup9ZO7ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1Zo/Bn+TIJOuSDCbZlOSSMWpOSfJyko3d7Y9H7HsmyaNd+8BUD0CSND5zeqjZAVxWVRuSHASsT3JvVW0eVfdAVZ25m2N8uKp+OqmeSpKmxB7P+KtqS1Vt6O6/AgwCC/d2xyRJe8e41viT9AHLgYfG2L0iyd8m+W6SpSPaC/h+kvVJ1ky4p5KkKdHLUg8ASQ4EbgUuraqto3ZvAN5VVb9McgZwO3B0t29lVT2f5O3AvUker6r7xzj+GmANwOLFiycwFElSL3o6408yl+HQv7mqbhu9v6q2VtUvu/t3A3OTHN5tP999fQH4S+DEsb5HVV1fVf1V1b9gwYIJDUaStGe9XNUT4AZgsKqu3k3Nb3d1JDmxO+5LSeZ3bwiTZD7wMeCxqeq8JGn8elnqWQmcCzyaZGPXdjmwGKCqrgXOAi5KsgP4FXB2VVWSI4C/7J4T5gD/raq+N8VjkCSNwx6Dv6oeBLKHmj8F/nSM9qeB4ybcO0nSlPOTu5LUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4Jakxewz+JEcmWZdkMMmmJJeMUXNKkpeTbOxufzxi3+okTyR5MslnpnoAkqTxmdNDzQ7gsqrakOQgYH2Se6tq86i6B6rqzJENSfYDvg58FBgCHk5y5xiPlSTtI3s846+qLVW1obv/CjAILOzx+CcCT1bV01X1D8AtwMcn2llJ0uT1csa/S5I+YDnw0Bi7VyT5W+B54N9V1SaGnyCeG1EzBHxgQj3twef+ahObn9+6tw4vSXvVe9/5Vq783aV7/fv0HPxJDgRuBS6tqtHpugF4V1X9MskZwO3A0UDGOFTt5vhrgDUAixcv7rVbkqRx6in4k8xlOPRvrqrbRu8f+URQVXcnuSbJ4Qyf4R85onQRw68IXqeqrgeuB+jv7x/zyWFP9sUzpSTNdL1c1RPgBmCwqq7eTc1vd3UkObE77kvAw8DRSZYkeQtwNnDnVHVekjR+vZzxrwTOBR5NsrFruxxYDFBV1wJnARcl2QH8Cji7qgrYkeQPgHuA/YAbu7V/SdI0yXA+v7n09/fXwMDAdHdDkmaMJOurqr+XWj+5K0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSY/YY/EmOTLIuyWCSTUkueYPaE5K8luSsEW2vJdnY3e6cqo5LkiZmTg81O4DLqmpDkoOA9UnurarNI4uS7Ad8Ebhn1ON/VVXLpqa7kqTJ2uMZf1VtqaoN3f1XgEFg4RilfwjcCrwwpT2UJE2pca3xJ+kDlgMPjWpfCPwecO0YD5uXZCDJ3yT5lxPspyRpivSy1ANAkgMZPqO/tKq2jtr9FeCPquq1JKMfuriqnk9yFPDXSR6tqqfGOP4aYA3A4sWLxzMGSdI49HTGn2Quw6F/c1XdNkZJP3BLkmeAs4Brdp7dV9Xz3dengR8w/Irhdarq+qrqr6r+BQsWjHcckqQe9XJVT4AbgMGqunqsmqpaUlV9VdUHfAf4t1V1e5JDkuzfHedwYCWweaxjSJL2jV6WelYC5wKPJtnYtV0OLAaoqrHW9Xf6HeC6JL9h+EnmC6OvBpIk7Vt7DP6qehB43cL9G9SfP+L+/wT+2YR6JknaK/zkriQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWrMHoM/yZFJ1iUZTLIpySVvUHtCkteSnDWi7bwkP+5u501VxyVJEzOnh5odwGVVtSHJQcD6JPdW1eaRRUn2A74I3DOi7VDgSqAfqO6xd1bVz6dsBJKkcdnjGX9VbamqDd39V4BBYOEYpX8I3Aq8MKLtdODeqvpZF/b3Aqsn3WtJ0oSNa40/SR+wHHhoVPtC4PeAa0c9ZCHw3IjtIcZ+0pAk7SM9B3+SAxk+o7+0qraO2v0V4I+q6rXRDxvjULWb469JMpBk4MUXX+y1W5KkcepljZ8kcxkO/Zur6rYxSvqBW5IAHA6ckWQHw2f4p4yoWwT8YKzvUVXXA9cD9Pf3j/nkIEmavD0Gf4bT/AZgsKquHqumqpaMqL8JuKuqbu/e3P3PSQ7pdn8M+Oykey1JmrBezvhXAucCjybZ2LVdDiwGqKrR6/q7VNXPkvwn4OGu6T9W1c8m0V9J0iTtMfir6kHGXqvfXf35o7ZvBG4cd88kSXuFn9yVpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMakqqa7D6+T5EXg7yf48MOBn05hd2YCxzz7tTZecMzj9a6qWtBL4Zsy+CcjyUBV9U93P/Ylxzz7tTZecMx7k0s9ktQYg1+SGjMbg//66e7ANHDMs19r4wXHvNfMujV+SdIbm41n/JKkNzBrgj/J6iRPJHkyyWemuz9TJcmRSdYlGUyyKcklXfuhSe5N8uPu6yFde5J8tfs5PJLk+OkdwcQl2S/JD5Pc1W0vSfJQN+ZvJXlL175/t/1kt79vOvs9UUneluQ7SR7v5nvFbJ/nJGu73+vHkvxFknmzbZ6T3JjkhSSPjWgb97wmOa+r/3GS8ybTp1kR/En2A74O/HPgvcA5Sd47vb2aMjuAy6rqd4APAp/uxvYZ4L6qOhq4r9uG4Z/B0d1tDfBn+77LU+YSYHDE9heBP+nG/HPgU137p4CfV9U/Bf6kq5uJ/gvwvao6BjiO4bHP2nlOshC4GOivqvcB+wFnM/vm+SZg9ai2cc1rkkOBK4EPACcCV+58spiQqprxN2AFcM+I7c8Cn53ufu2lsd4BfBR4AnhH1/YO4Inu/nXAOSPqd9XNpBuwqPsHcSpwFxCGP9gyZ/ScA/cAK7r7c7q6TPcYxjnetwJ/N7rfs3megYXAc8Ch3bzdBZw+G+cZ6AMem+i8AucA141o///qxnubFWf8/OMv0E5DXdus0r20XQ48BBxRVVsAuq9v78pmy8/iK8C/B37TbR8G/KKqdnTbI8e1a8zd/pe7+pnkKOBF4Jvd8tafJ5nPLJ7nqvo/wFXAs8AWhudtPbN7nnca77xO6XzPluDPGG2z6nKlJAcCtwKXVtXWNyodo21G/SySnAm8UFXrRzaPUVo97Jsp5gDHA39WVcuBV/nHl/9jmfFj7pYqPg4sAd4JzGd4qWO02TTPe7K7MU7p2GdL8A8BR47YXgQ8P019mXJJ5jIc+jdX1W1d80+SvKPb/w7gha59NvwsVgL/IskzwC0ML/d8BXhbkjldzchx7Rpzt/9g4Gf7ssNTYAgYqqqHuu3vMPxEMJvn+SPA31XVi1W1HbgN+BCze553Gu+8Tul8z5bgfxg4ursa4C0Mv0F05zT3aUokCXADMFhVV4/YdSew85398xhe+9/Z/vvd1QEfBF7e+ZJypqiqz1bVoqrqY3gu/7qq/jWwDjirKxs95p0/i7O6+hl1JlhV/xd4Lsl7uqbTgM3M4nlmeInng0l+q/s93znmWTvPI4x3Xu8BPpbkkO6V0se6tomZ7jc9pvDNkzOAHwFPAf9huvszheM6ieGXdI8AG7vbGQyvbd4H/Lj7emhXH4avcHoKeJThKyamfRyTGP8pwF3d/aOA/w08Cfx3YP+ufV63/WS3/6jp7vcEx7oMGOjm+nbgkNk+z8DngMeBx4D/Cuw/2+YZ+AuG38PYzvCZ+6cmMq/ABd3YnwT+zWT65Cd3Jakxs2WpR5LUI4Nfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TG/D9d/rgKtrw+dQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mask = np.array([1., 1., 0., 1/2])\n",
    "w_t = -2 * mask * P_tilde[[-1]] @ x_t_tilde\n",
    "plt.plot(w_t.ravel(), label='wage')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2.5, 2.5, 2.5, ..., 2.5, 2.5, 2.5]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w_t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4\n",
    "\n",
    "**Note:** \n",
    "- Changing the notation is not necessary because you're not reusing it for something else.\n",
    "- Your $\\tilde{R}$ in the pseudo-code is not symmetric\n",
    "- Shouldn't your $W$ matrix be $\\left[\\begin{array}{cc}\n",
    "\\boldsymbol{0} & 2q\\end{array}\\right]$ in the pseudo-code?\n",
    "- The `LQ` class does minimization so all the signs need to be reversed\n",
    "- Shouldn't $\\tilde{R}_f$ be symmetric?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 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}\n",
    "0 & 0 & 0 & 0\\\\\n",
    "0 & 0 & 0 & -\\frac{1}{2}\\\\\n",
    "0 & 0 & 0 & 0\\\\\n",
    "0 & -\\frac{1}{2} & 0 & f\n",
    "\\end{array}\\right]$ and $\\tilde{Q}_{f}=\\left[\\begin{array}{c}\n",
    "d\\end{array}\\right]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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:\n",
    "\n",
    "$$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)$$\n",
    "\n",
    "I'll also solve your Sylvester equation first in case I'm wrong. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "R_tilde_f = np.array([[0., 0., 0., 0.],\n",
    "                      [0., 0., 0., -1/2],\n",
    "                      [0., 0., 0., 0.],\n",
    "                      [0., -1/2, 0., f],])\n",
    "\n",
    "Q_tilde_f = np.array([[d]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SciPy has a solver for Sylvester equations of the form: \n",
    "\n",
    "$$AX+XB=Q$$ \n",
    "\n",
    "Therefore, we want to find matrices $A_{S}$, $B_{S}$ and $Q_{S}$ such that:\n",
    "\n",
    "$$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)$$ \n",
    "\n",
    "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: \n",
    "\n",
    "$$\\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}$$\n",
    "\n",
    "Post-multiplying by $\\left(A_{p}-B_{p}F_{p}\\right)$, we obtain the required equation.\n",
    "\n",
    "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)$. \n",
    "\n",
    "**Note:** The inverses don't necessarily exists depending on parameter values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "ename": "LinAlgError",
     "evalue": "Singular matrix",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mLinAlgError\u001b[0m                               Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-10-905271a33b0f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0mscipy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0msolve_sylvester\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mA_sylvester\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mβ\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mA_tilde\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      4\u001b[0m \u001b[0mB_sylvester\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mA_tilde\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mB_tilde\u001b[0m \u001b[0;34m@\u001b[0m \u001b[0mF_tilde\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mQ_sylvester\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mA_sylvester\u001b[0m \u001b[0;34m@\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mR_tilde_f\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mF_tilde\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m \u001b[0;34m@\u001b[0m \u001b[0mQ_tilde_f\u001b[0m \u001b[0;34m@\u001b[0m \u001b[0mF_tilde\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36minv\u001b[0;34m(a)\u001b[0m\n\u001b[1;32m    530\u001b[0m     \u001b[0msignature\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'D->D'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misComplexType\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m'd->d'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    531\u001b[0m     \u001b[0mextobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_linalg_error_extobj\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 532\u001b[0;31m     \u001b[0mainv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_umath_linalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msignature\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msignature\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mextobj\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mextobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    533\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mwrap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mainv\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult_t\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcopy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    534\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/anaconda3/lib/python3.7/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36m_raise_linalgerror_singular\u001b[0;34m(err, flag)\u001b[0m\n\u001b[1;32m     87\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     88\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 89\u001b[0;31m     \u001b[0;32mraise\u001b[0m \u001b[0mLinAlgError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Singular matrix\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     90\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     91\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_raise_linalgerror_nonposdef\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mLinAlgError\u001b[0m: Singular matrix"
     ]
    }
   ],
   "source": [
    "from scipy.linalg import solve_sylvester\n",
    "\n",
    "A_sylvester = -np.linalg.inv(β * A_tilde.T)\n",
    "B_sylvester = A_tilde - B_tilde @ F_tilde\n",
    "Q_sylvester = A_sylvester @ (R_tilde_f + F_tilde.T @ Q_tilde_f @ F_tilde)\n",
    "\n",
    "P_sylv = solve_sylvester(A_sylvester, B_sylvester, Q_sylvester)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "P_sylv = np.eye(4)\n",
    "max_iter = 1000\n",
    "δ = 1\n",
    "for i in range(max_iter):\n",
    "    δ = (R_tilde_f + F_tilde.T @ Q_tilde_f @ F_tilde + β * A_tilde.T @ P_sylv @ (A_tilde - B_tilde @ F_tilde)) - P_sylv\n",
    "    P_sylv += δ\n",
    "    \n",
    "    if np.max(np.abs(δ)) < 1e-12:\n",
    "        break\n",
    "        \n",
    "if i == max_iter - 1:\n",
    "    print('Max number of iterations reached.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[25., 25., 25., ..., 25., 25., 25.]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "-2 * P_sylv[[-1]] @ x_t_tilde"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[False, False, False, ..., False, False, False]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "-2 * P_sylv[[-1]] @ x_t_tilde == w_t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This doesn't recover the wage computed manually above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from quantecon import solve_discrete_lyapunov\n",
    "\n",
    "A_lyapunov = np.sqrt(β) * (A_tilde - B_tilde @ F_tilde)\n",
    "B_lyapunov = R_tilde_f + F_tilde.T @ Q_tilde_f @ F_tilde\n",
    "P_w_2 = solve_discrete_lyapunov(A_lyapunov, B_lyapunov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2.5, 2.5, 2.5, ..., 2.5, 2.5, 2.5]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "-2 * P_w_2[[-1]] @ x_t_tilde"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ True,  True,  True, ...,  True,  True,  True]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check against the one that I computed by hand\n",
    "-2 * P_w_2[[-1]] @ x_t_tilde == w_t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we look ready to try more complex parametrizations.\n",
    "\n",
    "**TODO:** Refactor the code to avoid repetition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 274,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "β = 0.9\n",
    "\n",
    "α_s = 5.\n",
    "ρ_1 = 0.\n",
    "σ_s = 0.\n",
    "\n",
    "α_b = 0.\n",
    "ρ_2 = 0.\n",
    "σ_b = 0.\n",
    "\n",
    "f = 1.\n",
    "h = 1.\n",
    "q = 0.\n",
    "d = 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 275,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Planner LQ problem\n",
    "R_tilde = np.array([[0., 0., 0., 0.],\n",
    "                    [0., 0., 0., -1/2],\n",
    "                    [0., 0., 0., 1/2],\n",
    "                    [0., -1/2, 1/2, f + h + 4 * q]])\n",
    "\n",
    "Q_tilde = np.array([[d + q]])\n",
    "\n",
    "N_tilde = np.array([[0., 0., 0., 2 * q]])\n",
    "\n",
    "A_tilde = np.array([[1., 0., 0., 0.],\n",
    "                    [α_s, ρ_1, 0., 0.],\n",
    "                    [α_b, 0., ρ_2, 0.],\n",
    "                    [0., 0., 0., 1.]])\n",
    "\n",
    "B_tilde = np.array([[0.],\n",
    "                    [0.],\n",
    "                    [0.],\n",
    "                    [1.]])\n",
    "\n",
    "C_tilde = np.array([[0., 0.],\n",
    "                    [σ_s, 0.],\n",
    "                    [0., σ_b],\n",
    "                    [0., 0.]])\n",
    "\n",
    "lq_tilde = qe.LQ(Q_tilde, R_tilde, A_tilde, B_tilde, C=C_tilde, N=N_tilde, beta=β)\n",
    "P_tilde, F_tilde, d_tilde = lq_tilde.stationary_values()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the following:\n",
    "\n",
    "$$\\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)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 276,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.220446049250313e-16"
      ]
     },
     "execution_count": 276,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "flow_utility = R_tilde + F_tilde.T @ Q_tilde @ F_tilde - 2 * F_tilde.T @ N_tilde\n",
    "continuation_val = (A_tilde - B_tilde @ F_tilde).T @ P_tilde @ (A_tilde - B_tilde @ F_tilde)\n",
    "\n",
    "# Should be very close to 0 if above holds\n",
    "np.max(np.abs(flow_utility + β * continuation_val - P_tilde))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 277,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.684341886080802e-14"
      ]
     },
     "execution_count": 277,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Should be very close to 0\n",
    "np.max(np.abs(solve_discrete_lyapunov(A=np.sqrt(β) * (A_tilde - B_tilde @ F_tilde).T,\n",
    "                               B=R_tilde + F_tilde.T @ Q_tilde @ F_tilde - 2 * F_tilde.T @ N_tilde)\n",
    "              - P_tilde))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 278,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Break down the problem for the firm\n",
    "R_tilde_f = np.array([[0., 0., 0., 0.],\n",
    "                      [0., 0., 0., -1/2],\n",
    "                      [0., 0., 0., 0.],\n",
    "                      [0., -1/2, 0., f],])\n",
    "\n",
    "Q_tilde_f = np.array([[d]])\n",
    "\n",
    "A_lyapunov = np.sqrt(β) * (A_tilde - B_tilde @ F_tilde).T\n",
    "B_lyapunov = R_tilde_f + F_tilde.T @ Q_tilde_f @ F_tilde\n",
    "P_w = solve_discrete_lyapunov(A=A_lyapunov, B=B_lyapunov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 279,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Same approach for the worker\n",
    "R_tilde_w = np.array([[0., 0., 0., 0.],\n",
    "                      [0., 0., 0., 0.],\n",
    "                      [0., 0., 0., 1/2],\n",
    "                      [0., 0., 1/2, h + 4 * q]])\n",
    "\n",
    "Q_tilde_w = np.array([[q]])\n",
    "\n",
    "N_tilde_w = np.array([[0., 0., 0., 2 * q]])\n",
    "\n",
    "B_lyapunov2 = R_tilde_w + F_tilde.T @ Q_tilde_w @ F_tilde - 2 * F_tilde.T @ N_tilde_w\n",
    "\n",
    "P_w2 = solve_discrete_lyapunov(A=A_lyapunov, B=B_lyapunov2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 280,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ True,  True,  True,  True],\n",
       "       [ True,  True,  True,  True],\n",
       "       [ True,  True,  True,  True],\n",
       "       [ True,  True,  True,  True]])"
      ]
     },
     "execution_count": 280,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R_tilde_w + R_tilde_f == R_tilde"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 281,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ True]])"
      ]
     },
     "execution_count": 281,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q_tilde_w + Q_tilde_f == Q_tilde"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 282,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ True,  True,  True,  True]])"
      ]
     },
     "execution_count": 282,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N_tilde_w == N_tilde"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 283,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.684341886080802e-14"
      ]
     },
     "execution_count": 283,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Should be very close to 0\n",
    "P_w3 = solve_discrete_lyapunov(A=A_lyapunov, B=B_lyapunov+B_lyapunov2)\n",
    "np.max(np.abs(P_w3 - P_tilde))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Implied wages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 284,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[3.3864]])"
      ]
     },
     "execution_count": 284,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "- 2 * P_w[[-1]] @ x_0_tilde.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 285,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[3.3864]])"
      ]
     },
     "execution_count": 285,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2 * P_w2[[-1]] @ x_0_tilde.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 286,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.]])"
      ]
     },
     "execution_count": 286,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2 * P_w3[[-1]] @ x_0_tilde.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do value functions add up?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 287,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-31.25]])"
      ]
     },
     "execution_count": 287,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_0_tilde @ P_w @ x_0_tilde.T + x_0_tilde @ P_w2 @ x_0_tilde.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 288,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-31.25]])"
      ]
     },
     "execution_count": 288,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_0_tilde @ P_w3 @ x_0_tilde.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 289,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.7729,  1.    ,  0.    , -3.4183]])"
      ]
     },
     "execution_count": 289,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Wage function with d=1 (computed by hand)\n",
    "(np.array([0., 1., 0., -2]) - 2 * F_tilde)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 290,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2.5]])"
      ]
     },
     "execution_count": 290,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(np.array([0., 1., 0., -2]) - 2 * F_tilde) @ x_0_tilde.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 291,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Firm LQ problem\n",
    "S_w = -2 * P_w[[-1]]\n",
    "S_s = np.array([[0., 1., 0., 0.]])\n",
    "cross_term = -(S_s - S_w) / 2\n",
    "\n",
    "n = A_tilde.shape[0]\n",
    "\n",
    "R_hat = np.block([[np.zeros((n, n)), cross_term.T],\n",
    "                  [cross_term, f]])\n",
    "\n",
    "Q_hat = np.array([[d]])\n",
    "\n",
    "N_hat = np.zeros((1, n + 1))\n",
    "\n",
    "A_hat = np.block([[A_tilde - B_tilde @ F_tilde, np.zeros((n, 1))],\n",
    "                  [np.zeros((1, n)), np.ones((1, 1))]])\n",
    "\n",
    "B_hat = np.array([[0.],\n",
    "                  [0.],\n",
    "                  [0.],\n",
    "                  [0.],\n",
    "                  [1.]])\n",
    "\n",
    "C_hat = np.vstack((C_tilde, np.zeros((1, C_tilde.shape[1]))))\n",
    "\n",
    "lq_hat = qe.LQ(Q_hat, R_hat, A_hat, B_hat, C=C_hat, N=N_hat, beta=β)\n",
    "P_hat, F_hat, d_hat = lq_hat.stationary_values()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 292,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.8864,  0.    ,  0.    ,  0.7092]])"
      ]
     },
     "execution_count": 292,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F_tilde"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 293,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.2292,  0.    ,  0.    , -0.1964,  0.5884]])"
      ]
     },
     "execution_count": 293,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F_hat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 294,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAF2JJREFUeJzt3X+QVeWd5/H3l6Zjr4BGARMjSEvGVUAQsMWgxpA1g2gyJChWYpkRxsz4I+NOJhkr0eyWOG4sswm7OtZKjGUcMmMWdYnZVaPRlcUyUTcIuvEHgqCiaU0EdfyBhgk/nv3j3osN9L2nu7nt5bm8X1Vd995zTp/zPH3w49Pfc/o5kVJCktRcBjS6AZKk+jPcJakJGe6S1IQMd0lqQoa7JDUhw12SmpDhLklNyHCXpCZkuEtSExrYqAMPGzYstbe3N+rwkpSlFStWvJZSGl60XcPCvb29neXLlzfq8JKUpYh4sSfbWZaRpCZkuEtSEzLcJakJGe6S1IQMd0lqQoa7JDUhw12SmlDD7nPvqwW33clBL93NvfvO5O2WDze6OZLUa2M/th/z/mxcvx4ju5H78E3rmL3xv7P/tjcb3RRJ2mNlN3I/85iR8DzMP/No+MjYRjdHkvZI2Y3c35ca3QBJ2mPlF+4RpddkuEtSNfmFO9HoBkjSHi/DcK9w5C5J1eQX7pZlJKlQfuEuSSqUYbhXau6O3CWpmvzC3bKMJBXKL9wlSYUyDHfLMpJUJL9wtywjSYXyC3dJUqEMw92yjCQVyS/cw+kHJKlIfuFe4cBdkqrKMNwty0hSkfzC3bKMJBUqDPeIuCki1kfEU1XWnx0RT5S/Ho6Io+vfzG54K6QkVdWTkftCYEaN9S8An0opTQD+E3BDHdpVg2UZSSpS+AzVlNKDEdFeY/3DXT7+X2DE7jerBqsyklSo3jX3rwD31Hmf3bMsI0lVFY7ceyoiPk0p3E+ssc15wHkAhx56aF+PVH413CWpmrqM3CNiAnAj8PmU0uvVtksp3ZBS6kgpdQwfPryvB+vb90nSXmS3wz0iDgVuB/48pfTs7jephyzLSFJVhWWZiFgETAOGRUQnMA9oBUgpXQ9cBgwFFkRpVL0lpdTRXw22LCNJxXpyt8xZBev/EvjLurWoiGUZSSqU31+oVliWkaSqMgx3yzKSVCS/cLcsI0mF8gv3CssyklRVhuFuWUaSiuQX7pZlJKlQfuFeYVlGkqrKMNwty0hSkfzC3bKMJBXKL9wrLMtIUlUZhrtlGUkqkl+4W5aRpEL5hXuFZRlJqirDcLcsI0lF8gt3yzKSVCi/cK9w4C5JVWUY7pZlJKlIfuFuWUaSCuUX7hXeLSNJVWUY7pZlJKlIfuFuVUaSCuUX7hWWZSSpqgzD3bKMJBXJL9y9W0aSCuUX7hWWZSSpqgzD3bKMJBXJL9wty0hSofzCvcKyjCRVlWG4W5aRpCL5hbtlGUkqlF+4V1iWkaSqMgx3yzKSVCS/cLcsI0mF8gv3CssyklRVhuFuWUaSihSGe0TcFBHrI+KpKuuPjIhHIuJfI+Li+jdxlwP2+yEkKXc9GbkvBGbUWP8G8DfA/Ho0qMcsy0hSVYXhnlJ6kFKAV1u/PqX0KLC5ng2rzrKMJBXJr+ZuWUaSCn2g4R4R50XE8ohYvmHDht3bmWUZSarqAw33lNINKaWOlFLH8OHD+7gXR+6SVMSyjCQ1oYFFG0TEImAaMCwiOoF5QCtASun6iPgosBzYD9gWEX8LjE0pvd1vrS4dvF93L0k5Kwz3lNJZBet/D4yoW4sKebeMJBWxLCNJTSi/cK+wLCNJVeUb7pZlJKmq/MLdsowkFcov3Cssy0hSVRmGu3fLSFKR/MLdsowkFcov3Cssy0hSVRmGu2UZSSqSX7hblpGkQvmFe4VlGUmqqnBumT2PZRmpyObNm+ns7GTTpk2Nbor6qK2tjREjRtDa2tqn788v3C3LSIU6OzsZMmQI7e3thP/NZCelxOuvv05nZyeHHXZYn/ZhWUZqQps2bWLo0KEGe6YigqFDh+7Wb14ZhrtlGaknDPa87e75yy/c/QcrZW3atGksX7680c0A4LLLLuP++++vy77a29t57bXXADj++ON79b2XX3458+fPr0s7KvKruVdYlpHUS1u3bqWlpWX75yuuuKJfjvPwww/3y357I7+Ru2UZaY+3bt06jjzySObMmcOECROYPXs277333i7bXXjhhXR0dDBu3DjmzZu3fXl7ezvz5s1j8uTJjB8/nlWrVgGlEe65557LtGnTGD16NNdee+3277n55puZMmUKEydO5Pzzz2fr1q0ADB48mMsuu4zjjjuORx55ZIfjz507l8WLF28/5re//W2mTp1KR0cHjz32GKeccgof//jHuf766wF44IEHOOmkk5g1axZjx47lggsuYNu2bbv0a/Dgwdvff//73+fYY49lwoQJO/Txyiuv5IgjjuAzn/kMq1ev7vXPuEh+I3fLMlKv/P2dT7Pylfo+0njsx/Zj3p+Nq7nN6tWr+dGPfsQJJ5zAueeey4IFC7j44ot32ObKK6/kwAMPZOvWrZx88sk88cQTTJgwAYBhw4bx2GOPsWDBAubPn8+NN94IwKpVq1i6dCnvvPMORxxxBBdeeCFr167l1ltv5aGHHqK1tZWvfvWr/OQnP+Gcc87h3Xff5aijjurRKH3kyJE88sgjfP3rX2fu3Lk89NBDbNq0iXHjxnHBBRcAsGzZMlauXMmoUaOYMWMGt99+O7Nnz+52f/fddx9r1qxh2bJlpJSYOXMmDz74IIMGDeKWW27h8ccfZ8uWLUyePJljjjmmsH29keHIvcyyjLRHGzlyJCeccAIAX/7yl/nVr361yza33XYbkydPZtKkSTz99NOsXLly+7rTTz8dgGOOOYZ169ZtX/7Zz36WffbZh2HDhnHQQQfx6quvsmTJElasWMGxxx7LxIkTWbJkCc8//zwALS0tnHHGGT1q88yZMwEYP348xx13HEOGDGH48OG0tbXx5ptvAjBlyhRGjx5NS0sLZ511Vrf9qrjvvvu47777mDRpEpMnT2bVqlWsWbOGX/7yl8yaNYt9992X/fbbb/tx6ym/kbtlGalXikbY/WXnuz12/vzCCy8wf/58Hn30UQ444ADmzp27w61/++yzD1AK5y1btuyyvOu6lBJz5szhqquu2qUdbW1tO9TZa6nse8CAATscZ8CAAdvbUNSvrlJKXHrppZx//vk7LL/mmmv6/W6m/EbulmWkLLz00kvba9yLFi3ixBNP3GH922+/zaBBg9h///159dVXueeee/p8rJNPPpnFixezfv16AN544w1efPHFvje+hmXLlvHCCy+wbds2br311l361dUpp5zCTTfdxMaNGwF4+eWXWb9+PSeddBI/+9nP+MMf/sA777zDnXfeWfd2ZjhyL7MsI+3RxowZw49//GPOP/98Dj/8cC688MId1h999NFMmjSJcePGMXr06O0lnL4YO3Ys3/nOd5g+fTrbtm2jtbWV6667jlGjRu1uN3YxdepULrnkEp588sntF1ermT59Os888wxTp04FShdab775ZiZPnswXv/hFJk6cyKhRo/jkJz9Z93ZGalBIdnR0pD7d67pxA8z/EzhtPkz5q/o3TGoCzzzzDGPGjGnY8detW8fnPvc5nnrqqYa1oT888MADzJ8/n7vuuusDOV535zEiVqSUOoq+17KMJDUhyzKS6q69vb3pRu1Q+uvaadOmNboZPZLfyN27ZSSpUH7hbllGkgrlF+4VlmUkqap8w92yjCRVlV+4W5aRsnDttdcyZswYzj77bO644w6++93v1m3fzTxVb714t4ykfrFgwQLuueee7Y+J627+lC1btjBwYO9jqJmn6q2X/Ebu3i0j7fEuuOACnn/+eWbOnMnVV1/NwoULueiii4DSNLvf+MY3+PSnP823vvUtLr/8cubMmcP06dNpb2/n9ttv55vf/Cbjx49nxowZbN68eZf9N/NUvfWS38jdsozUO/dcAr9/sr77/Oh4OLV6meX666/nF7/4BUuXLmXYsGEsXLhwh/XPPvss999/Py0tLVx++eU899xzLF26lJUrVzJ16lR++tOf8r3vfY9Zs2bx85//nC984Qs1m9NMU/XWS37hXmFZRsrWmWeeucNMjaeeeiqtra2MHz+erVu3MmPGDKA09W7X6X6r6TpV78aNGxkyZAhDhgzpdqpeYPtUvbXCvTJVL8DGjRtZs2YN77zzzvapersed0+UYbhblpF6pcYIu1EGDRq0w+euU+22trZunw6361S7tTTTVL31Ulhzj4ibImJ9RHT7t8RRcm1ErI2IJyJicv2bucMB+3X3kppTLlP11ktPRu4Lgf8G/FOV9acCh5e/jgN+UH7tX5ZlJPVCLlP11kuPpvyNiHbgrpTSUd2s+yHwQEppUfnzamBaSul3tfbZ5yl//3UjXHUI/OkVcMLXev/90l6g0VP+7mk+6Kl666XRU/4eAvy2y+fO8rL+YVlGkgrV44Jqd2nb7a8DEXEecB7AoYceuntHtSwjqYdymqq3Xuoxcu8ERnb5PAJ4pbsNU0o3pJQ6Ukodw4cP7+PhvFtGkorUI9zvAM4p3zXzCeCtonr7brEsI/VIox6hqfrY3fNXWJaJiEXANGBYRHQC84DW8sGvB+4GTgPWAu8Bf7FbLeop/+FKVbW1tfH6668zdOjQbO7L1vtSSrz++uu0tbX1eR+F4Z5SOqtgfQL+us8t6DXLMlKRESNG0NnZyYYNGxrdFPVRW1sbI0aM6PP35/cXqlGuJDlyl6pqbW3dPhuj9k75zQpZ+RXTcJekqjIM90qTDXdJqibfcE+7zsUsSSrJMNwrZRnDXZKqyS/cAQhr7pJUQ57hHgMcuUtSDZmGexjuklRDpuE+AO+WkaTq8g13R+6SVFWe4Y5lGUmqJc9wjwHeLSNJNWQa7t4KKUm1ZBruXlCVpFoyDXdr7pJUS57h7gVVSaopz3D3gqok1ZRxuDtyl6RqMg13yzKSVEum4e7dMpJUS57h7gVVSaopz3D3gqok1WS4S1ITyjTcLctIUi35hrsXVCWpqkzD3fvcJamWPMPdu2UkqaY8w90LqpJUU6bh7shdkmrJNNz9C1VJqiXfcHfkLklV5RnuXlCVpJryDHcvqEpSTYa7JDWhTMMdyzKSVEOm4e7dMpJUS4/CPSJmRMTqiFgbEZd0s35URCyJiCci4oGIGFH/pu5wREfuklRDYbhHRAtwHXAqMBY4KyLG7rTZfOCfUkoTgCuAq+rd0B0b5a2QklRLT0buU4C1KaXnU0p/BG4BPr/TNmOBJeX3S7tZX19eUJWkmnoS7ocAv+3yubO8rKvfAGeU388ChkTE0J13FBHnRcTyiFi+YcOGvrS3siNH7pJUQ0/CPbpZtvOw+WLgUxHxOPAp4GVgyy7flNINKaWOlFLH8OHDe93Y91vkBVVJqmVgD7bpBEZ2+TwCeKXrBimlV4DTASJiMHBGSumtejVyF9bcJammnozcHwUOj4jDIuJDwJeAO7puEBHDIqKyr0uBm+rbzJ2FNXdJqqEw3FNKW4CLgHuBZ4DbUkpPR8QVETGzvNk0YHVEPAt8BLiyn9pb4gVVSaqpJ2UZUkp3A3fvtOyyLu8XA4vr27QaImDb1g/scJKUm0z/QtW7ZSSplkzD3btlJKmWfMPdsowkVZVnuA8YCMlwl6Rq8g33bbv8jZQkqSzTcG+xLCNJNWQa7o7cJakWw12SmlCe4R4thrsk1ZBnuA8YaM1dkmrINNwduUtSLZmGuyN3Saol43B35C5J1WQc7o7cJamaTMPdmrsk1ZJpuFuWkaRaDHdJakL5hjsJtvnADknqTqbh3lJ6dfQuSd3KNNzLj3413CWpW4a7JDUhw12SmlCe4d7SWnrd+sfGtkOS9lB5hnvrvym9bv5DY9shSXuovMN9y6bGtkOS9lB5hvvAysj9vca2Q5L2UHmGu2UZSaop03Dft/S62bKMJHUn03BvK71alpGkbmUa7uWRuxdUJalbeYb7hwaXXje91dh2SNIeKs9w33coEPDua41uiSTtkfIM95aBsO+BsPHVRrdEkvZIeYY7wOCPwtsvN7oVkrRHyjfcD54AL6+ALc4vI0k7G9iTjSJiBvAPQAtwY0rpuzutPxT4MfDh8jaXpJTurnNbdzRmJvxmEfyXfwv7DoNo0P+nIhpzXEn5mvTncPxF/XqIwnCPiBbgOuBPgU7g0Yi4I6W0sstm/xG4LaX0g4gYC9wNtPdDe9935Gkw+x9h7RL448Z+PVR1qUHHlZS1wQf1+yF6MnKfAqxNKT0PEBG3AJ8HuoZ7AvYrv98feKWejazqqNNLX5KkHfQk3A8Bftvlcydw3E7bXA7cFxH/HhgEfKYurZMk9UlPCtXdFZV3rkecBSxMKY0ATgP+OWLXInhEnBcRyyNi+YYNG3rfWklSj/Qk3DuBkV0+j2DXsstXgNsAUkqPAG3AsJ13lFK6IaXUkVLqGD58eN9aLEkq1JNwfxQ4PCIOi4gPAV8C7thpm5eAkwEiYgylcHdoLkkNUhjuKaUtwEXAvcAzlO6KeToiroiImeXN/g74q4j4DbAImJtS8lYSSWqQHt3nXr5n/e6dll3W5f1K4IT6Nk2S1Ff5/oWqJKkqw12SmlA0qjQeERuAF/v47cOAvW2+X/u8d7DPe4fd6fOolFLh7YYNC/fdERHLU0odjW7HB8k+7x3s897hg+izZRlJakKGuyQ1oVzD/YZGN6AB7PPewT7vHfq9z1nW3CVJteU6cpck1ZBduEfEjIhYHRFrI+KSRrenXiJiZEQsjYhnIuLpiPhaefmBEfG/I2JN+fWA8vKIiGvLP4cnImJyY3vQNxHREhGPR8Rd5c+HRcSvy/29tTyfERGxT/nz2vL69ka2u68i4sMRsTgiVpXP9dS94Bx/vfxv+qmIWBQRbc12niPipohYHxFPdVnW6/MaEXPK26+JiDm706aswr3LU6FOBcYCZ5Wf/NQMtgB/l1IaA3wC+Oty3y4BlqSUDgeWlD9D6WdwePnrPOAHH3yT6+JrlOYsqvjPwNXl/v4LpRlHKb/+S0rpT4Cry9vl6B+AX6SUjgSOptT3pj3HEXEI8DdAR0rpKEqP4fwSzXeeFwIzdlrWq/MaEQcC8yg9L2MKMK/yP4Q+SSll8wVMBe7t8vlS4NJGt6uf+vq/KD3acDVwcHnZwcDq8vsfAmd12X77drl8UZo+egnw74C7KD074DVg4M7nm9LEdVPL7weWt4tG96GX/d0PeGHndjf5Oa487OfA8nm7CzilGc8zpUeLPtXX80rpuRg/7LJ8h+16+5XVyJ3unwp1SIPa0m/Kv4pOAn4NfCSl9DuA8mvl4YvN8LO4BvgmsK38eSjwZirNRAo79ml7f8vr3ypvn5PRlKbC/sdyKerGiBhEE5/jlNLLwHxK04L/jtJ5W0Fzn+eK3p7Xup7v3MK9J0+FylpEDAZ+CvxtSuntWpt2syybn0VEfA5Yn1Ja0XVxN5umHqzLxUBgMvCDlNIk4F3e/1W9O9n3uVxW+DxwGPAxSo/hPLWbTZvpPBep1se69j23cO/JU6GyFRGtlIL9Jyml28uLX42Ig8vrDwbWl5fn/rM4AZgZEeuAWyiVZq4BPhwRlamou/Zpe3/L6/cH3vggG1wHnUBnSunX5c+LKYV9s55jKD1P+YWU0oaU0mbgduB4mvs8V/T2vNb1fOcW7j15KlSWIiKAHwHPpJT+a5dVdwCVq+ZzKNXiK8vPKV95/wTwVuVXwByklC5NKY1IKbVTOo//J6V0NrAUmF3ebOf+Vn4Os8vbZzWiSyn9HvhtRBxRXnQysJImPcdlLwGfiIh9y//GK31u2vPcRW/P673A9Ig4oPwbz/Tysr5p9EWIPly0OA14FngO+A+Nbk8d+3UipV/BngD+X/nrNEr1xiXAmvLrgeXtg9KdQ88BT1K6G6Hh/ehj36cBd5XfjwaWAWuB/wHsU17eVv68trx+dKPb3ce+TgSWl8/z/wQOaPZzDPw9sAp4CvhnYJ9mO8+UnkD3O2AzpRH4V/pyXoFzy31fC/zF7rTJv1CVpCaUW1lGktQDhrskNSHDXZKakOEuSU3IcJekJmS4S1ITMtwlqQkZ7pLUhP4/tTFN9yzQQPgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_0_hat = np.array([[1., 5., 0., 1.25, 1.25]]).T\n",
    "x_t_hat = lq_hat.compute_sequence(x_0_hat, ts_length=1000)[0]\n",
    "\n",
    "plt.plot(x_t_hat[3], label='planner implied')\n",
    "plt.plot(x_t_hat[4], label='firm implied');\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 295,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.25, 1.25, 1.25, ..., 1.25, 1.25, 1.25])"
      ]
     },
     "execution_count": 295,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_t_hat[3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 296,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 296,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.var(x_t_hat[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 297,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.25  , 0.9892, 0.8819, ..., 0.8068, 0.8068, 0.8068])"
      ]
     },
     "execution_count": 297,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_t_hat[4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 298,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0002357102288177184"
      ]
     },
     "execution_count": 298,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.var(x_t_hat[4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 299,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.44322084543531903"
      ]
     },
     "execution_count": 299,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Sup norm\n",
    "np.max(np.abs(x_t_hat[3] - x_t_hat[4]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 300,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Worker LQ problem\n",
    "S_b = np.array([[0., 0., 1., 0.]])\n",
    "cross_term = (S_b - S_w) / 2\n",
    "\n",
    "n = A_tilde.shape[0]\n",
    "\n",
    "R_hat = np.block([[np.zeros((n, n)), cross_term.T],\n",
    "                  [cross_term, h + 4 * q]])\n",
    "\n",
    "Q_hat = np.array([[q]])\n",
    "\n",
    "N_hat = np.array([[0., 0., 0., 0., 2 * q]])\n",
    "\n",
    "A_hat = np.block([[A_tilde - B_tilde @ F_tilde, np.zeros((n, 1))],\n",
    "                  [np.zeros((1, n)), np.ones((1, 1))]])\n",
    "\n",
    "B_hat = np.array([[0.],\n",
    "                  [0.],\n",
    "                  [0.],\n",
    "                  [0.],\n",
    "                  [1.]])\n",
    "\n",
    "C_hat = np.vstack((C_tilde, np.zeros((1, C_tilde.shape[1]))))\n",
    "\n",
    "lq_hat = qe.LQ(Q_hat, R_hat, A_hat, B_hat, C=C_hat, N=N_hat, beta=β)\n",
    "P_hat, F_hat, d_hat = lq_hat.stationary_values()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 301,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGLVJREFUeJzt3X10VdWd//H3F4hkBGoRokOLJVCpPEiIIQgRgqEI0ipYWhylZQGlrYjTX2f82fGhv7WIMrqsFVvGaS1DKeKMbbRaSq0zWkYQKdSKiTKKPAgWdFIZebIVEDo8fH9/5OQ2QE5OHm5y2ZfPa60scs/ZOWefnKwPO9997o65OyIikl3aZboDIiKSfgp3EZEspHAXEclCCncRkSykcBcRyUIKdxGRLKRwFxHJQgp3EZEspHAXEclCHTJ14u7du3t+fn6mTi8iEqSqqqo97p6X1C5j4Z6fn09lZWWmTi8iEiQze7sx7VSWERHJQgp3EZEslBjuZrbYzHaZ2YaY/f9gZuujjw1mdszMzk1/V0VEpLEaM3JfAoyP2+nu97t7obsXAncAL7j7vjT1T0REmiEx3N19NdDYsJ4CVLSoRyIi0mJpq7mb2dnUjPB/3kCbG8ys0swqd+/ena5Ti4jISdI5oToBWNtQScbdF7p7sbsX5+UlPqYpIiLNlM7n3K+nrUoy638K+7a3yalERNLuE8PhwjGteoq0hLuZnQNcDkxNx/EadPR/Ydns2jO3+ulERNJu5N9nPtzNrAIoA7qbWTVQDuQAuPuCqNkkYLm7H2ylfv6FH6/5d0w5lP7fVj+diEiIEsPd3ac0os0Sah6ZbAPeNqcREQmY3qEqIpKFwg13U71dRCROeOHuKsuIiCQJL9xFRCRRgOFeO3JXWUZEJE6A4S4iIknCDXdNqIqIxAov3DWhKiKSKLxwFxGRRAGGuyZURUSSBBjuIiKSJNxw14SqiEis8MJdE6oiIonCC3cREUkUcLirLCMiEifAcFdZRkQkSYDhHtGEqohIrPDCXROqIiKJwgt3ERFJFHC4qywjIhInwHBXWUZEJEmA4R7RhKqISKzwwl0TqiIiicILdxERSRRwuKssIyISJ+BwFxGROOGGuyZURURihRfumlAVEUkUXriLiEiigMNdZRkRkTgBhrvKMiIiSQIM94gmVEVEYoUX7ppQFRFJFF64i4hIIoW7iEgWSgx3M1tsZrvMbEMDbcrMbL2ZvWFmL6S3iydTWUZEJEljRu5LgPFxO83so8BDwER3Hwhcm56uiYhIcyWGu7uvBvY10OSLwFJ3fydqvytNfYvrUM2/elpGRCRWOmrunwK6mtkqM6sys2lpOKaIiLRAhzQdYwgwBvgr4EUz+527v3lyQzO7AbgB4BOf+EQLT6uRu4hInHSM3KuBZ939oLvvAVYDg+tr6O4L3b3Y3Yvz8vKaeTpNqIqIJElHuP8SKDWzDmZ2NjAM2JSG44qISDMllmXMrAIoA7qbWTVQDuQAuPsCd99kZs8CrwHHgUXuHvvYZNpoQlVEJFZiuLv7lEa0uR+4Py09Sj5Zm5xGRCRkAb9DVSN3EZE4AYa7Ru4iIkkCDHcREUkSbrhrQlVEJFZ44a4JVRGRROGFe4pG7iIicQIMd43cRUSSBBjuIiKSJNxw14SqiEis8MJdE6oiIonCC/cUjdxFROIEGO4auYuIJAkw3EVEJEm44a4JVRGRWOGFuyZURUQShRfuKRq5i4jECTDcNXIXEUkSYLiLiEiScMNdE6oiIrHCC3dNqIqIJAov3FM0chcRiRNguGvkLiKSJMBwFxGRJOGGuyZURURihRfumlAVEUkUXrinaOQuIhIn4HAXEZE4CncRkSwUbrhrQlVEJFZ44a4JVRGRROGFu4iIJAo43FWWERGJE2C4qywjIpIkwHCPaEJVRCRWeOGuCVURkUThhbuIiCTqkNTAzBYDVwO73P3ievaXAb8Etkeblrr73HR2UkSa5siRI1RXV3P48OFMd0WaKTc3l549e5KTk9Osr08Md2AJ8H3gXxto8xt3v7pZPWgylWVEklRXV9OlSxfy8/MxzU8Fx93Zu3cv1dXV9O7du1nHSCzLuPtqYF+zjt6a9AMrEuvw4cN069ZNwR4oM6Nbt24t+s0rXTX3EjP7LzN7xswGxjUysxvMrNLMKnfv3t28M2lCVaRRFOxha+n9S0e4vwL0cvfBwD8Dy+IauvtCdy929+K8vLw0nFpEQlNWVkZlZWWmuwHAnDlzeO6559JyrPz8fPbs2QPAZZdd1qSvvfPOO5k3b15a+lGrMTX3Brn7B3U+/w8ze8jMurv7npYeu2EalYhI0xw7doz27dunXs+d2zrPfvz2t79tleM2RYtH7mb21xb9/mBml0bH3NvS48ZTWUbkdLdjxw769evH9OnTKSgoYPLkyXz44YentJs9ezbFxcUMHDiQ8vLy1Pb8/HzKy8spKipi0KBBbN68GagZ4c6cOZOysjL69OnDgw8+mPqaRx99lEsvvZTCwkJmzZrFsWPHAOjcuTNz5sxh2LBhvPjiiyecf8aMGTz55JOpc37rW9+ipKSE4uJiXnnlFa688ko++clPsmDBAgBWrVrFqFGjmDRpEgMGDODGG2/k+PHjp1xX586dU5/ff//9DB06lIKCghOu8Z577uGiiy7iiiuuYMuWLU3+HidpzKOQFUAZ0N3MqoFyIAfA3RcAk4HZZnYUOARc794GhXHVE0Ua5a5fvcHGdz9IbtgEAz72EconxE6vAbBlyxZ+/OMfM2LECGbOnMlDDz3EN7/5zRPa3HPPPZx77rkcO3aMMWPG8Nprr1FQUABA9+7deeWVV3jooYeYN28eixYtAmDz5s08//zz7N+/n4suuojZs2ezbds2Hn/8cdauXUtOTg433XQTP/nJT5g2bRoHDx7k4osvbtQo/YILLuDFF1/k5ptvZsaMGaxdu5bDhw8zcOBAbrzxRgDWrVvHxo0b6dWrF+PHj2fp0qVMnjy53uMtX76crVu3sm7dOtydiRMnsnr1ajp16sRjjz3Gq6++ytGjRykqKmLIkCGJ/WuKxHB39ykJ+79PzaOSbUMTqiJBuOCCCxgxYgQAU6dO5cEHHzwl3H/2s5+xcOFCjh49ys6dO9m4cWMq3D//+c8DMGTIEJYuXZr6mquuuoqOHTvSsWNHzjvvPN577z1WrFhBVVUVQ4cOBeDQoUOcd955ALRv354vfOELjerzxIkTARg0aBAHDhygS5cudOnShdzcXP74xz8CcOmll9KnTx8ApkyZwpo1axoM9+XLl3PJJZcAcODAAbZu3cr+/fuZNGkSZ5999gnnTacW19xF5PSWNMJuLSc/7XHy6+3btzNv3jxefvllunbtyowZM0549K9jx45ATTgfPXr0lO1197k706dP59577z2lH7m5uSfU2RtSe+x27dqdcJ527dql+pB0XXW5O3fccQezZs06Yfv8+fNb/WmmgJcfUFlG5HT2zjvvpGrcFRUVjBw58oT9H3zwAZ06deKcc87hvffe45lnnmn2ucaMGcOTTz7Jrl27ANi3bx9vv/128zvfgHXr1rF9+3aOHz/O448/fsp11XXllVeyePFiDhw4AMAf/vAHdu3axahRo/jFL37BoUOH2L9/P7/61a/S3s8AR+4qy4iEoH///jzyyCPMmjWLvn37Mnv27BP2Dx48mEsuuYSBAwfSp0+fVAmnOQYMGMDdd9/NuHHjOH78ODk5OfzgBz+gV69eLb2MU5SUlHD77bfz+uuvpyZX44wbN45NmzZRUlIC1Ey0PvrooxQVFXHddddRWFhIr169KC0tTXs/rS3mPutTXFzszXrW9b034IeXwd/8Kwy4Jv0dE8kCmzZton///hk7/44dO7j66qvZsGFDxvrQGlatWsW8efN4+umn2+R89d1HM6ty9+Kkrw2vLKMJVRGRRAGWZUTkdJefn591o3aoeXdtWVlZprvRKOGN3FM0oSoiEifAcFdZRkQkSYDhHtE7VEVEYoUX7ppQFRFJFF64i8gZpe5Suk21ZMkSvv71r5+y/amnnuLb3/52S7sGnLj42Fe/+lU2btzY6K9dtWoVV1/dOn/ELuCnZVSWEcl2tSs7ptvEiRNbZT2X2sXNTgcBjtxVlhE53X3nO99JLcd788038+lPfxqAFStWMHXqVKBmSYJBgwZx8cUXc9ttt6W+Nm6J3kOHDjF+/Hh+9KMfAfFL/D788MN86lOf4vLLL2ft2rX19q/uiH7GjBnMnj2b0aNH06dPH1544QVmzpxJ//79mTFjxgn9uuWWWygqKmLMmDHU99fk6v4hkuXLl1NSUkJRURHXXnttagmCZ599ln79+jFy5MgTFkRLt3BH7ppQFWmcZ26H/3k9vcf860HwmfiyxqhRo3jggQf4xje+QWVlJX/+8585cuQIa9asobS0lHfffZfbbruNqqoqunbtyrhx41i2bBmf+9zn6l2i98CBA1x//fVMmzaNadOmsWnTpnqX+B07dizl5eVUVVVxzjnnMHr06NSKjA15//33WblyJU899RQTJkxg7dq1LFq0iKFDh7J+/XoKCws5ePAgRUVFPPDAA8ydO5e77rqL73+//gVx9+zZw913381zzz1Hp06duO+++/jud7/Lrbfeyte+9jVWrlzJhRdeyHXXXdf0730jBThyF5HT3ZAhQ6iqqmL//v107NiRkpISKisr+c1vfkNpaSkvv/wyZWVl5OXl0aFDB770pS+xevVqoP4leq+55hq+/OUvM23aNIATlvgtLCxkxYoV/P73v+ell15KHfess85qdHhOmDABM2PQoEGcf/75DBo0iHbt2jFw4EB27NgB1KwMWXu8qVOnsmbNmtjj/e53v2Pjxo2MGDGCwsJCHnnkEd5++202b95M79696du3L2aW+i2mNYQ3ctfTMiJN08AIu7Xk5OSQn5/Pww8/zGWXXUZBQQHPP/88b731Fv379+fNN9+M/dr6lugdMWIEzzzzDF/84hcxs9glfpctW9aspXQbs9TvyZKW+h07diwVFRUnbF+/fn2b/eHygEfuKsuInM5GjRrFvHnzGDVqFKWlpSxYsIDCwkLMjGHDhvHCCy+wZ88ejh07RkVFBZdffnnssebOnUu3bt246aabgPglfocNG8aqVavYu3cvR44c4Yknnkjb9Rw/fjz1VMxPf/rTBpf6HT58OGvXrmXbtm0AfPjhh7z55pv069eP7du389ZbbwGcEv7pFGC4a+QuEoLS0lJ27txJSUkJ559/Prm5uamlbXv06MG9997L6NGjGTx4MEVFRVxzTcOrvM6fP5/Dhw9z6623nrDEb0FBAWPHjmXnzp306NGDO++8k5KSEq644gqKiorSdj2dOnXijTfeYMiQIaxcuZI5c+bEts3Ly2PJkiVMmTKFgoIChg8fzubNm8nNzWXhwoVcddVVjBw5slWWJK4V3pK/774KC8tgymNw0WfS3i+RbJDpJX+zUefOnVNPvLSVM2vJXxERSRReuGtCVUQyoK1H7S0VXrinaEJVRCROgOGukbtIY2RqPk3So6X3L8BwF5Ekubm57N27VwEfKHdn79695ObmNvsY4b2JqZaWHxCJ1bNnT6qrq+td/0TCkJubS8+ePZv99eGFuwYiIolycnLo3bt3prshGRRwWUYjdxGROAGGu4buIiJJAgx3ERFJEm64a0JVRCRWeOGuR7tERBKFF+4pGrmLiMQJMNw1chcRSRJguIuISJJww11VGRGRWInhbmaLzWyXmW1IaDfUzI6Z2eT0da8emlAVEUnUmJH7EmB8Qw3MrD1wH/DrNPSpkTR0FxGJkxju7r4a2JfQ7P8APwd2paNTCT1q/VOIiASuxTV3M/s4MAlY0PLuiIhIOqRjQnU+cJu7H0tqaGY3mFmlmVW2eClSvUNVRCRWOpb8LQYes5qw7Q581syOuvuykxu6+0JgIUBxcXHz6iuaUBURSdTicHf31KLRZrYEeLq+YE8/jdxFROIkhruZVQBlQHczqwbKgRwAd1edXUTkNJQY7u4+pbEHc/cZLepN487S+qcQEQlcwO9QVVlGRCROeOGuCVURkUThhXuKRu4iInECDncREYkTYLirLCMikiTAcI9oQlVEJFZ44a4JVRGRROGFe4pG7iIicQIOdxERiRNguKssIyKSJMBwj2hCVUQkVnjhrglVEZFE4YV7ikbuIiJxAg53ERGJE2C4qywjIpIkwHCPaEJVRCRWeOGuCVURkUThhbuIiCQKONxVlhERiRNguKssIyKSJMBwj2hCVUQkVnjhrglVEZFE4YW7iIgkCjjcVZYREYkTYLirLCMikiTAcI9oQlVEJFa44S4iIrHCC3dVZUREEoUX7ikqy4iIxAkw3DV0FxFJEmC4RzShKiISK9xwFxGRWOGFu5YfEBFJFF64p6gsIyISJ8Bw18hdRCRJgOEe0cBdRCRWYrib2WIz22VmG2L2X2Nmr5nZejOrNLOR6e+miIg0RWNG7kuA8Q3sXwEMdvdCYCawKA39iqcJVRGRRInh7u6rgX0N7D/gnkrcTrRZUVx1GRGROGmpuZvZJDPbDPw7NaP3VqSRu4hIkrSEu7v/wt37AZ8D/jGunZndENXlK3fv3t2yk+odqiIisdL6tExUwvmkmXWP2b/Q3YvdvTgvLy+dpxYRkTpaHO5mdqFZzTDazIqAs4C9LT1uLE2oiogk6pDUwMwqgDKgu5lVA+VADoC7LwC+AEwzsyPAIeC6OhOsrUhlGRGROInh7u5TEvbfB9yXth4l0shdRCRJwO9Q1chdRCROuOEuIiKxwgt3TaiKiCQKL9xTVJYREYkTYLhr5C4ikiTAcI9oQlVEJFa44S4iIrHCC3dNqIqIJAov3FNUlhERiRNwuIuISJwAw11lGRGRJAGGe0RPy4iIxAov3DWhKiKSKLxwT9HIXUQkTsDhLiIicQIMd5VlRESSJP6xjtPNE1XVXAvc8sR/UZ3zp0x3R0SkyQZ87COUTxjYqucIbuT+QU4eL+aWcqhdp0x3RUTktBXcyP0rU64Hrqck0x0RETmNBTdyFxGRZAp3EZEspHAXEclCCncRkSykcBcRyUIKdxGRLKRwFxHJQgp3EZEsZJ6hJXTNbDfwdjO/vDuwJ43dCYGu+cygaz4ztOSae7l7XlKjjIV7S5hZpbsXZ7ofbUnXfGbQNZ8Z2uKaVZYREclCCncRkSwUargvzHQHMkDXfGbQNZ8ZWv2ag6y5i4hIw0IduYuISAOCC3czG29mW8xsm5ndnun+pIuZXWBmz5vZJjN7w8z+Ltp+rpn9p5ltjf7tGm03M3sw+j68ZmZFmb2C5jGz9mb2qpk9Hb3ubWYvRdf7uJmdFW3vGL3eFu3Pz2S/m8vMPmpmT5rZ5uhel5wB9/jm6Gd6g5lVmFlutt1nM1tsZrvMbEOdbU2+r2Y2PWq/1cymt6RPQYW7mbUHfgB8BhgATDGzAZntVdocBW5x9/7AcOBvo2u7HVjh7n2BFdFrqPke9I0+bgB+2PZdTou/AzbVeX0f8L3oet8HvhJt/wrwvrtfCHwvaheifwKedfd+wGBqrj1r77GZfRz4BlDs7hcD7YHryb77vAQYf9K2Jt1XMzsXKAeGAZcC5bX/ITSLuwfzAZQAv67z+g7gjkz3q5Wu9ZfAWGAL0CPa1gPYEn3+L8CUOu1T7UL5AHpGP/SfBp4GjJo3dnQ4+X4DvwZKos87RO0s09fQxOv9CLD95H5n+T3+OPDfwLnRfXsauDIb7zOQD2xo7n0FpgD/Umf7Ce2a+hHUyJ2//KDUqo62ZZXoV9FLgJeA8919J0D073lRs2z4XswHbgWOR6+7AX9096PR67rXlLreaP+fovYh6QPsBh6OSlGLzKwTWXyP3f0PwDzgHWAnNfetiuy+z7Wael/Ter9DC3erZ1tWPe5jZp2BnwN/7+4fNNS0nm3BfC/M7Gpgl7tX1d1cT1NvxL5QdACKgB+6+yXAQf7yq3p9gr/mqKxwDdAb+BjQiZqyxMmy6T4nibvGtF57aOFeDVxQ53VP4N0M9SXtzCyHmmD/ibsvjTa/Z2Y9ov09gF3R9tC/FyOAiWa2A3iMmtLMfOCjZlb7h9vrXlPqeqP95wD72rLDaVANVLv7S9HrJ6kJ+2y9xwBXANvdfbe7HwGWApeR3fe5VlPva1rvd2jh/jLQN5ppP4uaiZmnMtyntDAzA34MbHL379bZ9RRQO2s+nZpafO32adHM+3DgT7W/AobA3e9w957unk/NfVzp7l8CngcmR81Ovt7a78PkqH1QIzp3/x/gv83somjTGGAjWXqPI+8Aw83s7OhnvPaas/Y+19HU+/prYJyZdY1+4xkXbWueTE9CNGPS4rPAm8BbwP/LdH/SeF0jqfkV7DVgffTxWWrqjSuArdG/50btjZonh94CXqfmaYSMX0czr70MeDr6vA+wDtgGPAF0jLbnRq+3Rfv7ZLrfzbzWQqAyus/LgK7Zfo+Bu4DNwAbg34CO2XafgQpq5hSOUDMC/0pz7iswM7r2bcCXW9InvUNVRCQLhVaWERGRRlC4i4hkIYW7iEgWUriLiGQhhbuISBZSuIuIZCGFu4hIFlK4i4hkof8PFK9BktHcM6oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_0_hat = np.array([[1., 5., 0., 1.25, 1.25]]).T\n",
    "x_t_hat = lq_hat.compute_sequence(x_0_hat, ts_length=1000)[0]\n",
    "\n",
    "plt.plot(x_t_hat[3], label='planner implied')\n",
    "plt.plot(x_t_hat[4], label='worked implied')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 302,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.44322084543531926"
      ]
     },
     "execution_count": 302,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Sup norm\n",
    "np.max(np.abs(x_t_hat[3] - x_t_hat[4]))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}