perm_income_cons.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='perm-income-cons'></a>\n",
    "<div id=\"qe-notebook-header\" align=\"right\" style=\"text-align:right;\">\n",
    "        <a href=\"https://quantecon.org/\" title=\"quantecon.org\">\n",
    "                <img style=\"width:250px;display:inline;\" width=\"250px\" src=\"https://assets.quantecon.org/img/qe-menubar-logo.svg\" alt=\"QuantEcon\">\n",
    "        </a>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimal Savings II: LQ Techniques\n",
    "\n",
    "\n",
    "<a id='index-1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Optimal Savings II: LQ Techniques](#Optimal-Savings-II:-LQ-Techniques)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [Setup](#Setup)  \n",
    "  - [The LQ Approach](#The-LQ-Approach)  \n",
    "  - [Implementation](#Implementation)  \n",
    "  - [Two Example Economies](#Two-Example-Economies)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Co-author: Chase Coleman**\n",
    "\n",
    "In addition to what’s in Anaconda, this lecture will need the following libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": true
   },
   "outputs": [],
   "source": [
    "!pip install --upgrade quantecon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "This lecture continues our analysis of the linear-quadratic (LQ) permanent income model of savings and consumption.\n",
    "\n",
    "As we saw in our [previous lecture](https://lectures.quantecon.org/perm_income.html) on this topic, Robert Hall [[Hal78]](https://lectures.quantecon.org/zreferences.html#hall1978) used the LQ permanent income model to restrict and interpret intertemporal comovements of nondurable consumption, nonfinancial income, and financial wealth.\n",
    "\n",
    "For example, we saw how the model asserts that for any covariance stationary process for nonfinancial income\n",
    "\n",
    "- consumption is a random walk  \n",
    "- financial wealth has a unit root and is cointegrated with consumption  \n",
    "\n",
    "\n",
    "Other applications use the same LQ framework.\n",
    "\n",
    "For example, a model isomorphic to the LQ permanent income model has been used by Robert Barro [[Bar79]](https://lectures.quantecon.org/zreferences.html#barro1979) to interpret intertemporal comovements of a government’s tax collections, its  expenditures net of debt service, and its public debt.\n",
    "\n",
    "This isomorphism means that in analyzing the LQ permanent income model, we are in effect also analyzing  the Barro tax smoothing model.\n",
    "\n",
    "It is just a matter of appropriately relabeling the variables in Hall’s model.\n",
    "\n",
    "In this lecture, we’ll\n",
    "\n",
    "- show how the solution to the LQ permanent income model can be obtained using LQ control methods.  \n",
    "- represent the model as a linear state space system as in [this lecture](https://lectures.quantecon.org/linear_models.html).  \n",
    "- apply [QuantEcon](http://quantecon.org/python_index.html)’s [LinearStateSpace](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lss.py) class to characterize statistical features of the consumer’s optimal consumption and borrowing plans.  \n",
    "\n",
    "\n",
    "We’ll then use these characterizations to construct a simple model of cross-section wealth and\n",
    "consumption dynamics in the spirit of Truman Bewley [[Bew86]](https://lectures.quantecon.org/zreferences.html#bewley86).\n",
    "\n",
    "(Later we’ll study other Bewley models—see [this lecture](https://lectures.quantecon.org/aiyagari.html))\n",
    "\n",
    "The model will prove useful for illustrating concepts such as\n",
    "\n",
    "- stationarity  \n",
    "- ergodicity  \n",
    "- ensemble moments and cross-section observations  \n",
    "\n",
    "\n",
    "Let’s start with some imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import quantecon as qe\n",
    "import numpy as np\n",
    "import scipy.linalg as la\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "Let’s recall the basic features of the model discussed in the [permanent income model](https://lectures.quantecon.org/perm_income.html).\n",
    "\n",
    "Consumer preferences are ordered  by\n",
    "\n",
    "\n",
    "<a id='equation-old1'></a>\n",
    "$$\n",
    "E_0 \\sum_{t=0}^\\infty \\beta^t u(c_t) \\tag{1}\n",
    "$$\n",
    "\n",
    "where $ u(c) = -(c - \\gamma)^2 $.\n",
    "\n",
    "The consumer maximizes [(1)](#equation-old1) by choosing a\n",
    "consumption, borrowing plan $ \\{c_t, b_{t+1}\\}_{t=0}^\\infty $\n",
    "subject to the sequence of budget constraints\n",
    "\n",
    "\n",
    "<a id='equation-old2'></a>\n",
    "$$\n",
    "c_t + b_t = \\frac{1}{1 + r} b_{t+1}  + y_t,\n",
    "\\quad t \\geq 0 \\tag{2}\n",
    "$$\n",
    "\n",
    "and the no-Ponzi condition\n",
    "\n",
    "\n",
    "<a id='equation-old42'></a>\n",
    "$$\n",
    "E_0 \\sum_{t=0}^\\infty \\beta^t b_t^2 < \\infty \\tag{3}\n",
    "$$\n",
    "\n",
    "The interpretation of all variables and parameters are the same as in the\n",
    "[previous lecture](https://lectures.quantecon.org/perm_income.html).\n",
    "\n",
    "We continue to assume that $ (1 + r) \\beta = 1 $.\n",
    "\n",
    "The dynamics of $ \\{y_t\\} $ again follow the linear state space model\n",
    "\n",
    "\n",
    "<a id='equation-sprob15ab2'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    z_{t+1} & = A z_t + C w_{t+1}\n",
    "    \\\\\n",
    "    y_t & = U  z_t\n",
    "\\end{aligned} \\tag{4}\n",
    "$$\n",
    "\n",
    "The restrictions on the shock process and parameters are the same as in our [previous lecture](https://lectures.quantecon.org/perm_income.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Digression on a Useful Isomorphism\n",
    "\n",
    "The LQ permanent income model of consumption is mathematically isomorphic with a version of\n",
    "Barro’s [[Bar79]](https://lectures.quantecon.org/zreferences.html#barro1979) model of tax smoothing.\n",
    "\n",
    "In the LQ permanent income model\n",
    "\n",
    "- the household faces an exogenous process of nonfinancial income  \n",
    "- the household wants to smooth consumption across states and time  \n",
    "\n",
    "\n",
    "In the Barro tax smoothing model\n",
    "\n",
    "- a government faces an exogenous sequence of government purchases (net of  interest payments on its debt)  \n",
    "- a government wants to smooth tax collections across states and time  \n",
    "\n",
    "\n",
    "If we set\n",
    "\n",
    "- $ T_t $, total tax collections in Barro’s model to consumption $ c_t $ in the LQ permanent income model.  \n",
    "- $ G_t $, exogenous government expenditures in Barro’s model to nonfinancial income $ y_t $ in the permanent income model.  \n",
    "- $ B_t $, government risk-free one-period assets falling due in Barro’s model to  risk-free one-period consumer debt $ b_t $ falling due in the LQ permanent income model.  \n",
    "- $ R $, the gross rate of return on risk-free one-period government debt in Barro’s model to the gross rate of return $ 1+r $ on financial assets in  the permanent income model of consumption.  \n",
    "\n",
    "\n",
    "then the two models are mathematically equivalent.\n",
    "\n",
    "All characterizations of a $ \\{c_t, y_t, b_t\\} $ in the LQ permanent income model automatically apply to a $ \\{T_t, G_t, B_t\\} $ process in the Barro model of tax smoothing.\n",
    "\n",
    "See [consumption and tax smoothing models](https://lectures.quantecon.org/smoothing.html) for further exploitation of an isomorphism between consumption and tax smoothing models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A Specification of the Nonfinancial Income Process\n",
    "\n",
    "For the purposes of this lecture, let’s assume $ \\{y_t\\} $ is a second-order univariate autoregressive process:\n",
    "\n",
    "$$\n",
    "y_{t+1} = \\alpha + \\rho_1 y_t + \\rho_2 y_{t-1} + \\sigma w_{t+1}\n",
    "$$\n",
    "\n",
    "We can map this into the linear state space framework in [(4)](#equation-sprob15ab2), as\n",
    "discussed in our lecture on  [linear models](https://lectures.quantecon.org/linear_models.html).\n",
    "\n",
    "To do so we take\n",
    "\n",
    "$$\n",
    "z_t =\n",
    "\\begin{bmatrix}\n",
    "    1 \\\\\n",
    "    y_t \\\\\n",
    "    y_{t-1}\n",
    "\\end{bmatrix},\n",
    "\\quad\n",
    "A = \\begin{bmatrix}\n",
    "          1 & 0 & 0 \\\\\n",
    "          \\alpha & \\rho_1 & \\rho_2  \\\\\n",
    "          0 & 1 & 0\n",
    "    \\end{bmatrix},\n",
    "\\quad\n",
    "C= \\begin{bmatrix}\n",
    "    0 \\\\\n",
    "    \\sigma \\\\\n",
    "    0\n",
    "    \\end{bmatrix},\n",
    "\\quad \\text{and} \\quad\n",
    "U = \\begin{bmatrix} 0 & 1 & 0 \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The LQ Approach\n",
    "\n",
    "[Previously](https://lectures.quantecon.org/perm_income.html#odr-pi) we solved the permanent income model  by solving a system of linear expectational difference equations subject to two boundary conditions.\n",
    "\n",
    "Here we solve the same model using [LQ methods](https://lectures.quantecon.org/lqcontrol.html) based on dynamic programming.\n",
    "\n",
    "After confirming that answers produced by the two methods agree, we apply [QuantEcon](http://quantecon.org/python_index.html)’s [LinearStateSpace](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lss.py)\n",
    "class to illustrate features of the model.\n",
    "\n",
    "Why solve a model in two distinct ways?\n",
    "\n",
    "Because by doing so we gather insights about the structure of the model.\n",
    "\n",
    "Our earlier approach based on solving a system of expectational difference equations brought to the fore the role of the consumer’s expectations about future nonfinancial income.\n",
    "\n",
    "On the other hand, formulating the model in terms of an LQ dynamic programming problem reminds us that\n",
    "\n",
    "- finding the state (of a dynamic programming problem) is an art, and  \n",
    "- iterations on a Bellman equation  implicitly jointly solve both a  forecasting problem and a control problem  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The LQ Problem\n",
    "\n",
    "Recall from our [lecture on LQ theory](https://lectures.quantecon.org/lqcontrol.html) that the optimal linear regulator problem is to choose\n",
    "a decision rule for $ u_t $ to minimize\n",
    "\n",
    "$$\n",
    "\\mathbb E\n",
    "\\sum_{t=0}^\\infty \\beta^t \\{x'_t R x_t+ u'_t Q u_t\\},\n",
    "$$\n",
    "\n",
    "subject to $ x_0 $ given and the law of motion\n",
    "\n",
    "\n",
    "<a id='equation-pilqsd'></a>\n",
    "$$\n",
    "x_{t+1} =  \\tilde A x_t+ \\tilde B u_t+ \\tilde C w_{t+1},\n",
    "\\qquad t\\geq 0, \\tag{5}\n",
    "$$\n",
    "\n",
    "where $ w_{t+1} $ is IID with mean vector zero and $ \\mathbb E w_t w'_t= I $.\n",
    "\n",
    "The tildes in $ \\tilde A, \\tilde B, \\tilde C $ are to avoid clashing with notation in [(4)](#equation-sprob15ab2).\n",
    "\n",
    "The value function for this problem is $ v(x) = - x'Px - d $, where\n",
    "\n",
    "- $ P $ is the unique positive semidefinite solution of the [corresponding matrix Riccati equation](https://lectures.quantecon.org/lqcontrol.html#riccati-equation).  \n",
    "- The scalar $ d $ is given by $ d=\\beta (1-\\beta)^{-1} {\\rm trace} ( P \\tilde C \\tilde C') $.  \n",
    "\n",
    "\n",
    "The optimal policy is $ u_t = -Fx_t $, where $ F := \\beta (Q+\\beta \\tilde B'P \\tilde B)^{-1} \\tilde B'P \\tilde A $.\n",
    "\n",
    "Under an optimal decision rule $ F $, the state vector $ x_t $ evolves according to $ x_{t+1} = (\\tilde A-\\tilde BF) x_t + \\tilde C w_{t+1} $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mapping into the LQ Framework\n",
    "\n",
    "To map into the LQ framework, we’ll use\n",
    "\n",
    "$$\n",
    "x_t :=\n",
    "    \\begin{bmatrix}\n",
    "        z_t \\\\\n",
    "        b_t\n",
    "    \\end{bmatrix} =\n",
    "    \\begin{bmatrix}\n",
    "        1 \\\\\n",
    "        y_t \\\\\n",
    "        y_{t-1} \\\\\n",
    "        b_t\n",
    "    \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "as the state vector and $ u_t := c_t - \\gamma $ as the control.\n",
    "\n",
    "With this notation and $ U_\\gamma := \\begin{bmatrix} \\gamma & 0 & 0\n",
    "\\end{bmatrix} $, we can write the state dynamics as in [(5)](#equation-pilqsd) when\n",
    "\n",
    "$$\n",
    "\\tilde A :=\n",
    "     \\begin{bmatrix}\n",
    "        A  &  0 \\\\\n",
    "        (1 + r)(U_\\gamma - U) & 1 + r\n",
    "     \\end{bmatrix}\n",
    "\\quad\n",
    "\\tilde B :=\n",
    "   \\begin{bmatrix}\n",
    "   0 \\\\\n",
    "   1 + r\n",
    "   \\end{bmatrix}\n",
    "\\quad \\text{and} \\quad\n",
    "\\tilde C :=\n",
    "   \\begin{bmatrix}\n",
    "   C \\\\ 0\n",
    "   \\end{bmatrix}\n",
    "   w_{t+1}\n",
    "$$\n",
    "\n",
    "Please confirm for yourself that, with these definitions, the LQ dynamics [(5)](#equation-pilqsd)  match the dynamics of $ z_t $ and $ b_t $ described above.\n",
    "\n",
    "To map utility into the quadratic form $ x_t' R x_t + u_t'Q u_t $ we can set\n",
    "\n",
    "- $ Q := 1 $ (remember that we are minimizing) and  \n",
    "- $ R := $ a $ 4 \\times 4 $ matrix of zeros  \n",
    "\n",
    "\n",
    "However, there is one problem remaining.\n",
    "\n",
    "We have no direct way to capture the non-recursive restriction [(3)](#equation-old42)\n",
    "on the debt sequence $ \\{b_t\\} $ from within the LQ framework.\n",
    "\n",
    "To try to enforce it, we’re going to use a trick: put a small penalty on $ b_t^2 $ in the criterion function.\n",
    "\n",
    "In the present setting, this means adding a small entry $ \\epsilon > 0 $ in the $ (4,4) $ position of $ R $.\n",
    "\n",
    "That will induce a (hopefully) small approximation error in the decision rule.\n",
    "\n",
    "We’ll check whether it really is small numerically soon."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "Let’s write some code to solve the model.\n",
    "\n",
    "One comment before we start is that the bliss level of consumption $ \\gamma $ in the utility function has no effect on the optimal decision rule.\n",
    "\n",
    "We saw this in the previous lecture  [permanent income](https://lectures.quantecon.org/perm_income.html).\n",
    "\n",
    "The reason is that it drops out of the Euler equation for consumption.\n",
    "\n",
    "In what follows we set it equal to unity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Exogenous Nonfinancial Income Process\n",
    "\n",
    "First, we create the objects for the optimal linear regulator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Set parameters\n",
    "α, β, ρ1, ρ2, σ = 10.0, 0.95, 0.9, 0.0, 1.0\n",
    "\n",
    "R = 1 / β\n",
    "A = np.array([[1., 0., 0.],\n",
    "              [α,  ρ1, ρ2],\n",
    "              [0., 1., 0.]])\n",
    "C = np.array([[0.], [σ], [0.]])\n",
    "G = np.array([[0., 1., 0.]])\n",
    "\n",
    "# Form LinearStateSpace system and pull off steady state moments\n",
    "μ_z0 = np.array([[1.0], [0.0], [0.0]])\n",
    "Σ_z0 = np.zeros((3, 3))\n",
    "Lz = qe.LinearStateSpace(A, C, G, mu_0=μ_z0, Sigma_0=Σ_z0)\n",
    "μ_z, μ_y, Σ_z, Σ_y = Lz.stationary_distributions()\n",
    "\n",
    "# Mean vector of state for the savings problem\n",
    "mxo = np.vstack([μ_z, 0.0])\n",
    "\n",
    "# Create stationary covariance matrix of x -- start everyone off at b=0\n",
    "a1 = np.zeros((3, 1))\n",
    "aa = np.hstack([Σ_z, a1])\n",
    "bb = np.zeros((1, 4))\n",
    "sxo = np.vstack([aa, bb])\n",
    "\n",
    "# These choices will initialize the state vector of an individual at zero\n",
    "# debt and the ergodic distribution of the endowment process. Use these to\n",
    "# create the Bewley economy.\n",
    "mxbewley = mxo\n",
    "sxbewley = sxo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to create the matrices for the LQ system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "A12 = np.zeros((3,1))\n",
    "ALQ_l = np.hstack([A, A12])\n",
    "ALQ_r = np.array([[0, -R, 0, R]])\n",
    "ALQ = np.vstack([ALQ_l, ALQ_r])\n",
    "\n",
    "RLQ = np.array([[0., 0., 0., 0.],\n",
    "                [0., 0., 0., 0.],\n",
    "                [0., 0., 0., 0.],\n",
    "                [0., 0., 0., 1e-9]])\n",
    "\n",
    "QLQ = np.array([1.0])\n",
    "BLQ = np.array([0., 0., 0., R]).reshape(4,1)\n",
    "CLQ = np.array([0., σ, 0., 0.]).reshape(4,1)\n",
    "β_LQ = β"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s print these out and have a look at them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "print(f\"A = \\n {ALQ}\")\n",
    "print(f\"B = \\n {BLQ}\")\n",
    "print(f\"R = \\n {RLQ}\")\n",
    "print(f\"Q = \\n {QLQ}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now create the appropriate instance of an LQ model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "lqpi = qe.LQ(QLQ, RLQ, ALQ, BLQ, C=CLQ, beta=β_LQ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We’ll save the implied optimal policy function soon compare them with what we get by\n",
    "employing an alternative solution method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "P, F, d = lqpi.stationary_values()  # Compute value function and decision rule\n",
    "ABF = ALQ - BLQ @ F  #  Form closed loop system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison with the Difference Equation Approach\n",
    "\n",
    "In our [first lecture](https://lectures.quantecon.org/perm_income.html) on the infinite horizon permanent\n",
    "income problem we used a different solution method.\n",
    "\n",
    "The method was based around\n",
    "\n",
    "- deducing the Euler equations that are the first-order conditions with respect to consumption and savings.  \n",
    "- using the budget constraints and boundary condition to complete a system of expectational linear difference equations.  \n",
    "- solving those equations to obtain the solution.  \n",
    "\n",
    "\n",
    "Expressed in state space notation, the solution took the form\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    z_{t+1} & = A z_t + C w_{t+1} \\\\\n",
    "    b_{t+1} & = b_t + U [ (I -\\beta A)^{-1} (A - I) ] z_t \\\\\n",
    "        y_t & = U z_t \\\\\n",
    "        c_t & = (1-\\beta) [ U (I-\\beta A)^{-1} z_t - b_t ]\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Now we’ll apply the formulas in this system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Use the above formulas to create the optimal policies for b_{t+1} and c_t\n",
    "b_pol = G @ la.inv(np.eye(3, 3) - β * A) @ (A - np.eye(3, 3))\n",
    "c_pol = (1 - β) * G @ la.inv(np.eye(3, 3) - β * A)\n",
    "\n",
    "# Create the A matrix for a LinearStateSpace instance\n",
    "A_LSS1 = np.vstack([A, b_pol])\n",
    "A_LSS2 = np.eye(4, 1, -3)\n",
    "A_LSS = np.hstack([A_LSS1, A_LSS2])\n",
    "\n",
    "# Create the C matrix for LSS methods\n",
    "C_LSS = np.vstack([C, np.zeros(1)])\n",
    "\n",
    "# Create the G matrix for LSS methods\n",
    "G_LSS1 = np.vstack([G, c_pol])\n",
    "G_LSS2 = np.vstack([np.zeros(1), -(1 - β)])\n",
    "G_LSS = np.hstack([G_LSS1, G_LSS2])\n",
    "\n",
    "# Use the following values to start everyone off at b=0, initial incomes zero\n",
    "μ_0 = np.array([1., 0., 0., 0.])\n",
    "Σ_0 = np.zeros((4, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`A_LSS` calculated as we have here should equal `ABF` calculated above using the LQ model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ABF - A_LSS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now compare pertinent elements of `c_pol` and `F`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "print(c_pol, \"\\n\", -F)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have verified that the two methods give the same solution.\n",
    "\n",
    "Now let’s create instances of the [LinearStateSpace](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lss.py) class and use it to do some interesting experiments.\n",
    "\n",
    "To do this, we’ll use the outcomes from our second method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two Example Economies\n",
    "\n",
    "In the spirit of Bewley models [[Bew86]](https://lectures.quantecon.org/zreferences.html#bewley86), we’ll generate panels of consumers.\n",
    "\n",
    "The examples differ only in  the initial states with which we endow the consumers.\n",
    "\n",
    "All other parameter values are kept the same in the two examples\n",
    "\n",
    "- In the first example, all consumers begin with zero nonfinancial income and zero debt.  \n",
    "  \n",
    "  - The consumers are thus *ex-ante* identical.  \n",
    "  \n",
    "- In the second example, while all begin with zero debt, we draw their initial income levels from the invariant distribution of financial income.  \n",
    "  \n",
    "  - Consumers are *ex-ante* heterogeneous.  \n",
    "  \n",
    "\n",
    "\n",
    "In the first example, consumers’ nonfinancial income paths  display\n",
    "pronounced transients early in the sample\n",
    "\n",
    "- these will affect outcomes in striking ways  \n",
    "\n",
    "\n",
    "Those transient effects will not be present in the second example.\n",
    "\n",
    "We use methods affiliated with the [LinearStateSpace](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lss.py) class to simulate the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First Set of Initial Conditions\n",
    "\n",
    "We generate  25 paths of the exogenous non-financial income process and the associated optimal consumption and debt paths.\n",
    "\n",
    "In the first set of graphs,  darker lines depict a particular sample path, while the lighter lines describe 24 other  paths.\n",
    "\n",
    "A second graph  plots a collection of simulations against the population distribution that we extract from the `LinearStateSpace` instance `LSS`.\n",
    "\n",
    "Comparing sample paths with population distributions at each date $ t $ is a useful exercise—see [our discussion](https://lectures.quantecon.org/lln_clt.html#lln-mr) of the laws of large numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "lss = qe.LinearStateSpace(A_LSS, C_LSS, G_LSS, mu_0=μ_0, Sigma_0=Σ_0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Population and Sample Panels\n",
    "\n",
    "In the code below, we use the [LinearStateSpace](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lss.py) class to\n",
    "\n",
    "- compute and plot population quantiles of the distributions of\n",
    "  consumption and debt for a population of consumers.  \n",
    "- simulate a group of 25 consumers and plot sample paths on the same\n",
    "  graph as the population distribution.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def income_consumption_debt_series(A, C, G, μ_0, Σ_0, T=150, npaths=25):\n",
    "    \"\"\"\n",
    "    This function takes initial conditions (μ_0, Σ_0) and uses the\n",
    "    LinearStateSpace class from QuantEcon to  simulate an economy\n",
    "    npaths times for T periods. It then uses that information to\n",
    "    generate some graphs related to the discussion below.\n",
    "    \"\"\"\n",
    "    lss = qe.LinearStateSpace(A, C, G, mu_0=μ_0, Sigma_0=Σ_0)\n",
    "\n",
    "    # Simulation/Moment Parameters\n",
    "    moment_generator = lss.moment_sequence()\n",
    "\n",
    "    # Simulate various paths\n",
    "    bsim = np.empty((npaths, T))\n",
    "    csim = np.empty((npaths, T))\n",
    "    ysim = np.empty((npaths, T))\n",
    "\n",
    "    for i in range(npaths):\n",
    "        sims = lss.simulate(T)\n",
    "        bsim[i, :] = sims[0][-1, :]\n",
    "        csim[i, :] = sims[1][1, :]\n",
    "        ysim[i, :] = sims[1][0, :]\n",
    "\n",
    "    # Get the moments\n",
    "    cons_mean = np.empty(T)\n",
    "    cons_var = np.empty(T)\n",
    "    debt_mean = np.empty(T)\n",
    "    debt_var = np.empty(T)\n",
    "    for t in range(T):\n",
    "        μ_x, μ_y, Σ_x, Σ_y = next(moment_generator)\n",
    "        cons_mean[t], cons_var[t] = μ_y[1], Σ_y[1, 1]\n",
    "        debt_mean[t], debt_var[t] = μ_x[3], Σ_x[3, 3]\n",
    "\n",
    "    return bsim, csim, ysim, cons_mean, cons_var, debt_mean, debt_var\n",
    "\n",
    "def consumption_income_debt_figure(bsim, csim, ysim):\n",
    "\n",
    "    # Get T\n",
    "    T =  bsim.shape[1]\n",
    "\n",
    "    # Create the first figure\n",
    "    fig, ax = plt.subplots(2, 1, figsize=(10, 8))\n",
    "    xvals = np.arange(T)\n",
    "\n",
    "    # Plot consumption and income\n",
    "    ax[0].plot(csim[0, :], label=\"c\", color=\"b\")\n",
    "    ax[0].plot(ysim[0, :], label=\"y\", color=\"g\")\n",
    "    ax[0].plot(csim.T, alpha=.1, color=\"b\")\n",
    "    ax[0].plot(ysim.T, alpha=.1, color=\"g\")\n",
    "    ax[0].legend(loc=4)\n",
    "    ax[0].set(title=\"Nonfinancial Income, Consumption, and Debt\",\n",
    "              xlabel=\"t\", ylabel=\"y and c\")\n",
    "\n",
    "    # Plot debt\n",
    "    ax[1].plot(bsim[0, :], label=\"b\", color=\"r\")\n",
    "    ax[1].plot(bsim.T, alpha=.1, color=\"r\")\n",
    "    ax[1].legend(loc=4)\n",
    "    ax[1].set(xlabel=\"t\", ylabel=\"debt\")\n",
    "\n",
    "    fig.tight_layout()\n",
    "    return fig\n",
    "\n",
    "def consumption_debt_fanchart(csim, cons_mean, cons_var,\n",
    "                              bsim, debt_mean, debt_var):\n",
    "    # Get T\n",
    "    T =  bsim.shape[1]\n",
    "\n",
    "    # Create percentiles of cross-section distributions\n",
    "    cmean = np.mean(cons_mean)\n",
    "    c90 = 1.65 * np.sqrt(cons_var)\n",
    "    c95 = 1.96 * np.sqrt(cons_var)\n",
    "    c_perc_95p, c_perc_95m = cons_mean + c95, cons_mean - c95\n",
    "    c_perc_90p, c_perc_90m = cons_mean + c90, cons_mean - c90\n",
    "\n",
    "    # Create percentiles of cross-section distributions\n",
    "    dmean = np.mean(debt_mean)\n",
    "    d90 = 1.65 * np.sqrt(debt_var)\n",
    "    d95 = 1.96 * np.sqrt(debt_var)\n",
    "    d_perc_95p, d_perc_95m = debt_mean + d95, debt_mean - d95\n",
    "    d_perc_90p, d_perc_90m = debt_mean + d90, debt_mean - d90\n",
    "\n",
    "\n",
    "    # Create second figure\n",
    "    fig, ax = plt.subplots(2, 1, figsize=(10, 8))\n",
    "    xvals = np.arange(T)\n",
    "\n",
    "    # Consumption fan\n",
    "    ax[0].plot(xvals, cons_mean, color=\"k\")\n",
    "    ax[0].plot(csim.T, color=\"k\", alpha=.25)\n",
    "    ax[0].fill_between(xvals, c_perc_95m, c_perc_95p, alpha=.25, color=\"b\")\n",
    "    ax[0].fill_between(xvals, c_perc_90m, c_perc_90p, alpha=.25, color=\"r\")\n",
    "    ax[0].set(title=\"Consumption/Debt over time\",\n",
    "              ylim=(cmean-15, cmean+15), ylabel=\"consumption\")\n",
    "\n",
    "    # Debt fan\n",
    "    ax[1].plot(xvals, debt_mean, color=\"k\")\n",
    "    ax[1].plot(bsim.T, color=\"k\", alpha=.25)\n",
    "    ax[1].fill_between(xvals, d_perc_95m, d_perc_95p, alpha=.25, color=\"b\")\n",
    "    ax[1].fill_between(xvals, d_perc_90m, d_perc_90p, alpha=.25, color=\"r\")\n",
    "    ax[1].set(xlabel=\"t\", ylabel=\"debt\")\n",
    "\n",
    "    fig.tight_layout()\n",
    "    return fig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let’s create figures with initial conditions of zero for $ y_0 $ and $ b_0 $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "out = income_consumption_debt_series(A_LSS, C_LSS, G_LSS, μ_0, Σ_0)\n",
    "bsim0, csim0, ysim0 = out[:3]\n",
    "cons_mean0, cons_var0, debt_mean0, debt_var0 = out[3:]\n",
    "\n",
    "consumption_income_debt_figure(bsim0, csim0, ysim0)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "consumption_debt_fanchart(csim0, cons_mean0, cons_var0,\n",
    "                          bsim0, debt_mean0, debt_var0)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is what is going on in the above graphs.\n",
    "\n",
    "For our simulation, we have set initial conditions $ b_0 = y_{-1} = y_{-2} = 0 $.\n",
    "\n",
    "Because $ y_{-1} = y_{-2} = 0 $, nonfinancial income $ y_t $ starts far below its stationary mean $ \\mu_{y, \\infty} $ and rises early in each simulation.\n",
    "\n",
    "Recall from  the [previous lecture](https://lectures.quantecon.org/perm_income.html) that we can represent the optimal decision rule for consumption in terms of the **co-integrating relationship**\n",
    "\n",
    "\n",
    "<a id='equation-old12'></a>\n",
    "$$\n",
    "(1-\\beta) b_t + c_t = (1-\\beta) E_t \\sum_{j=0}^\\infty \\beta^j y_{t+j} \\tag{6}\n",
    "$$\n",
    "\n",
    "So at time $ 0 $ we have\n",
    "\n",
    "$$\n",
    "c_0 = (1-\\beta) E_0 \\sum_{t=0}^\\infty \\beta^j y_{t}\n",
    "$$\n",
    "\n",
    "This tells us that consumption starts at the income that would be paid by an annuity whose value equals the expected discounted value of nonfinancial income at time $ t=0 $.\n",
    "\n",
    "To support that level of consumption, the consumer borrows a lot early and consequently builds up substantial debt.\n",
    "\n",
    "In fact, he or she incurs so much debt that eventually, in the stochastic steady state, he consumes less each period than his nonfinancial income.\n",
    "\n",
    "He uses the gap between consumption and nonfinancial income mostly to service the interest payments due on his debt.\n",
    "\n",
    "Thus, when we look at the panel of debt in the accompanying graph, we see that this is a group of *ex-ante* identical people each of whom starts with zero debt.\n",
    "\n",
    "All of them accumulate debt in anticipation of rising nonfinancial income.\n",
    "\n",
    "They expect their nonfinancial income to rise toward the invariant distribution of income, a consequence of our having started them at $ y_{-1} = y_{-2} = 0 $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cointegration Residual\n",
    "\n",
    "The following figure plots realizations of the left side of [(6)](#equation-old12), which,\n",
    "[as discussed in our last lecture](https://lectures.quantecon.org/perm_income.html#coint-pi), is called the **cointegrating residual**.\n",
    "\n",
    "As mentioned above, the right side can be thought of as an\n",
    "annuity payment on the expected present value of future income\n",
    "$ E_t \\sum_{j=0}^\\infty \\beta^j y_{t+j} $.\n",
    "\n",
    "Early along a realization, $ c_t $ is approximately constant while\n",
    "$ (1-\\beta) b_t $ and\n",
    "$ (1-\\beta) E_t \\sum_{j=0}^\\infty \\beta^j y_{t+j} $ both rise\n",
    "markedly as the household’s present value of income and borrowing rise\n",
    "pretty much together.\n",
    "\n",
    "This example illustrates the following point: the definition\n",
    "of cointegration implies that the cointegrating residual is\n",
    "*asymptotically* covariance stationary, not *covariance stationary*.\n",
    "\n",
    "The cointegrating residual for the specification with zero income and zero\n",
    "debt initially has a notable transient component that dominates its\n",
    "behavior early in the sample.\n",
    "\n",
    "By altering initial conditions, we shall remove this transient in our second example to be presented below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def cointegration_figure(bsim, csim):\n",
    "    \"\"\"\n",
    "    Plots the cointegration\n",
    "    \"\"\"\n",
    "    # Create figure\n",
    "    fig, ax = plt.subplots(figsize=(10, 8))\n",
    "    ax.plot((1 - β) * bsim[0, :] + csim[0, :], color=\"k\")\n",
    "    ax.plot((1 - β) * bsim.T + csim.T, color=\"k\", alpha=.1)\n",
    "\n",
    "    ax.set(title=\"Cointegration of Assets and Consumption\", xlabel=\"t\")\n",
    "\n",
    "    return fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "cointegration_figure(bsim0, csim0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A “Borrowers and Lenders” Closed Economy\n",
    "\n",
    "When we set $ y_{-1} = y_{-2} = 0 $ and $ b_0 =0 $ in the\n",
    "preceding exercise, we make debt “head north” early in the sample.\n",
    "\n",
    "Average debt in the cross-section rises and approaches the asymptote.\n",
    "\n",
    "We can regard these as outcomes of a “small open economy” that\n",
    "borrows from abroad at the fixed gross interest rate $ R = r+1 $ in\n",
    "anticipation of rising incomes.\n",
    "\n",
    "So with the economic primitives set as above, the economy converges to a\n",
    "steady state in which there is an excess aggregate supply of risk-free\n",
    "loans at a gross interest rate of $ R $.\n",
    "\n",
    "This excess supply is filled by “foreigner lenders” willing to make those loans.\n",
    "\n",
    "We can use virtually the same code to rig a “poor man’s Bewley [[Bew86]](https://lectures.quantecon.org/zreferences.html#bewley86) model” in the following way\n",
    "\n",
    "- as before, we start everyone at $ b_0 = 0 $.  \n",
    "- But instead of starting everyone at $ y_{-1} = y_{-2} = 0 $, we\n",
    "  draw $ \\begin{bmatrix} y_{-1} \\\\ y_{-2}   \\end{bmatrix} $ from\n",
    "  the invariant distribution of the $ \\{y_t\\} $ process.  \n",
    "\n",
    "\n",
    "This rigs a closed economy in which people are borrowing and lending\n",
    "with each other at a gross risk-free interest rate of\n",
    "$ R = \\beta^{-1} $.\n",
    "\n",
    "Across the group of people being analyzed, risk-free loans are in zero excess supply.\n",
    "\n",
    "We have arranged primitives so that $ R = \\beta^{-1} $ clears the market for risk-free loans at zero aggregate excess supply.\n",
    "\n",
    "So the risk-free loans are being made from one person to another within our closed set of agent.\n",
    "\n",
    "There is no need for foreigners to lend to our group.\n",
    "\n",
    "Let’s have a look at the corresponding figures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "out = income_consumption_debt_series(A_LSS, C_LSS, G_LSS, mxbewley, sxbewley)\n",
    "bsimb, csimb, ysimb = out[:3]\n",
    "cons_meanb, cons_varb, debt_meanb, debt_varb = out[3:]\n",
    "\n",
    "consumption_income_debt_figure(bsimb, csimb, ysimb)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "consumption_debt_fanchart(csimb, cons_meanb, cons_varb,\n",
    "                          bsimb, debt_meanb, debt_varb)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The graphs confirm the following outcomes:\n",
    "\n",
    "- As before, the consumption distribution spreads out over time.  \n",
    "\n",
    "\n",
    "But now there is some initial dispersion because there is *ex-ante* heterogeneity in the initial draws of $ \\begin{bmatrix} y_{-1} \\\\ y_{-2}   \\end{bmatrix} $.\n",
    "\n",
    "- As before, the cross-section distribution of debt spreads out over time.  \n",
    "- Unlike before, the average level of debt stays at zero, confirming that this is a closed borrower-and-lender economy.  \n",
    "- Now the cointegrating residual seems stationary, and not just asymptotically stationary.  \n",
    "\n",
    "\n",
    "Let’s have a look at the cointegration figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "cointegration_figure(bsimb, csimb)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "filename": "perm_income_cons.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Optimal Savings II: LQ Techniques"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}