re_with_feedback.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='re-with-feedback'></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": [
    "# Stability in Linear Rational Expectations Models\n",
    "\n",
    "\n",
    "<a id='index-0'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Stability in Linear Rational Expectations Models](#Stability-in-Linear-Rational-Expectations-Models)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [Linear difference equations](#Linear-difference-equations)  \n",
    "  - [Illustration: Cagan’s Model](#Illustration:-Cagan’s-Model)  \n",
    "  - [Some Python code](#Some-Python-code)  \n",
    "  - [Alternative code](#Alternative-code)  \n",
    "  - [Another perspective](#Another-perspective)  \n",
    "  - [Log money supply feeds back on log price level](#Log-money-supply-feeds-back-on-log-price-level)  \n",
    "  - [Big $ P $, little $ p $ interpretation](#Big-$-P-$,-little-$-p-$-interpretation)  \n",
    "  - [Fun with Sympy code](#Fun-with-Sympy-code)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to what’s in Anaconda, this lecture deploys the following libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": true
   },
   "outputs": [],
   "source": [
    "!pip install --upgrade quantecon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import quantecon as qe\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from sympy import *\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "This lecture studies stability in the context of an elementary rational expectations model.\n",
    "\n",
    "We study a rational expectations version of Philip Cagan’s model [[Cag56]](https://python-programming.quantecon.org/zreferences.html#cagan) linking\n",
    "the price level to the money supply.\n",
    "\n",
    "Cagan did not use a rational expectations version of his model, but Sargent [[Sar77]](https://python-programming.quantecon.org/zreferences.html#sargent77hyper) did.\n",
    "\n",
    "We study this model because it is intrinsically interesting and also because it  has a mathematical structure that\n",
    "also appears in virtually all  linear rational expectations model, namely, that a key  endogenous variable equals\n",
    "a mathematical expectation of a geometric sum of future values of another variable.\n",
    "\n",
    "In a rational expectations version of Cagan’s model, the endogenous variable is the price level or rate of inflation and\n",
    "the other variable is the money supply or the rate of change in the money supply.\n",
    "\n",
    "In this lecture, we’ll encounter:\n",
    "\n",
    "- a convenient formula for the expectation of geometric sum of future values of a variable  \n",
    "- a way of solving an expectational difference equation by mapping it into a vector first-order difference equation and appropriately manipulating an eigen decomposition of the transition matrix in order to impose stability  \n",
    "- a way to use a Big $ K $, little $ k $ argument to allow apparent feedback from endogenous to exogenous variables within a rational expectations equilibrium  \n",
    "- a use of eigenvector decompositions of matrices that allowed Blanchard and Khan (1981) [[BK80]](https://python-programming.quantecon.org/zreferences.html#blanchard-khan) and Whiteman (1983) [[Whi83]](https://python-programming.quantecon.org/zreferences.html#whiteman) to solve a class of linear rational expectations models  \n",
    "\n",
    "\n",
    "Cagan’s model with rational expectations\n",
    "is formulated as an **expectational difference equation** whose solution is a rational expectations equilibrium.\n",
    "\n",
    "We’ll start this lecture with a quick review of deterministic (i.e., non-random)\n",
    "first-order and second-order linear difference equations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear difference equations\n",
    "\n",
    "In this quick review of linear difference equations, we’ll use the *backward shift* or *lag* operator $ L $\n",
    "\n",
    "The lag operator $ L $  maps a sequence $ \\{x_t\\}_{t=0}^\\infty $ into the sequence $ \\{x_{t-1}\\}_{t=0}^\\infty $\n",
    "\n",
    "We’ll can use  $ L $ in linear difference equations by using the equality\n",
    "$ L x_t \\equiv x_{t-1} $ in algebraic expressions.\n",
    "\n",
    "Further,  the inverse $ L^{-1} $ of the lag operator is  the *forward shift*\n",
    "operator.\n",
    "\n",
    "In linear difference equations, we’ll often use the equaltiy  $ L^{-1} x_t \\equiv x_{t+1} $ in the the algebra\n",
    "below.\n",
    "\n",
    "The algebra of lag and forward shift operators often simplifies formulas for linear difference equations and their\n",
    "solutions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First order\n",
    "\n",
    "We want to solve a linear first-order scalar difference equation.\n",
    "\n",
    "First, let $ |\\lambda | < 1 $, and let\n",
    "$ \\{u_t\\}_{t=-\\infty}^\\infty $ be a bounded sequence of scalar real\n",
    "numbers.\n",
    "\n",
    "Let $ L $ be the lag operator defined by\n",
    "$ L x_t \\equiv x_{t-1} $ and let $ L^{-1} $ be the forward shift\n",
    "operator defined by $ L^{-1} x_t \\equiv x_{t+1} $.\n",
    "\n",
    "Then\n",
    "\n",
    "\n",
    "<a id='equation-equn-1'></a>\n",
    "$$\n",
    "(1 - \\lambda L) y_t = u_t, \\forall t \\tag{1}\n",
    "$$\n",
    "\n",
    "has solutions\n",
    "\n",
    "\n",
    "<a id='equation-equn-2'></a>\n",
    "$$\n",
    "y_t = (1 -\\lambda L)^{-1} u_t +k \\lambda^t \\tag{2}\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "$$\n",
    "y_t =  \\sum_{j=0}^\\infty \\lambda^j u_{t-j} +k \\lambda^t\n",
    "$$\n",
    "\n",
    "for any real number $ k $.\n",
    "\n",
    "You can verify this fact by applying $ (1-\\lambda L) $ to both sides\n",
    "of equation [(2)](#equation-equn-2) and noting that $ (1 - \\lambda L) \\lambda^t =0 $.\n",
    "\n",
    "To pin down $ k $ we need one condition imposed from outside (e.g.,\n",
    "an initial or terminal condition) on the path of $ y $.\n",
    "\n",
    "Now let $ | \\lambda | > 1 $.\n",
    "\n",
    "Rewrite equation [(1)](#equation-equn-1) as\n",
    "\n",
    "\n",
    "<a id='equation-equn-3'></a>\n",
    "$$\n",
    "y_{t-1} = \\lambda^{-1} y_t - \\lambda^{-1} u_t , \\forall t \\tag{3}\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "\n",
    "<a id='equation-equn-4'></a>\n",
    "$$\n",
    "(1 - \\lambda^{-1} L^{-1}) y_t = - \\lambda^{-1} u_{t+1}. \\tag{4}\n",
    "$$\n",
    "\n",
    "A solution is\n",
    "\n",
    "\n",
    "<a id='equation-equn-5'></a>\n",
    "$$\n",
    "y_t = - \\lambda^{-1}\\left({ 1 \\over  1 - \\lambda^{-1} L^{-1}} \\right)\n",
    "         u_{t+1} + k \\lambda^t \\tag{5}\n",
    "$$\n",
    "\n",
    "for any $ k $.\n",
    "\n",
    "To verify that this is a solution, check the consequences of operating\n",
    "on both sides of equation [(5)](#equation-equn-5) by $ (1 -\\lambda L) $ and compare to\n",
    "equation [(1)](#equation-equn-1).\n",
    "\n",
    "Solution [(2)](#equation-equn-2) exists for $ |\\lambda | < 1 $ because\n",
    "the distributed lag in $ u $ converges.\n",
    "\n",
    "Solution [(5)](#equation-equn-5) exists when $ |\\lambda| > 1 $ because the **distributed\n",
    "lead** in $ u $ converges.\n",
    "\n",
    "When $ |\\lambda | > 1 $, the distributed lag in $ u $ in [(2)](#equation-equn-2) may\n",
    "diverge, so that a solution of this form does not exist.\n",
    "\n",
    "The distributed lead in $ u $ in [(5)](#equation-equn-5) need not\n",
    "converge when $ |\\lambda| < 1 $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Second order\n",
    "\n",
    "Now consider the second order difference equation\n",
    "\n",
    "\n",
    "<a id='equation-equn-6'></a>\n",
    "$$\n",
    "(1-\\lambda_1 L) (1 - \\lambda_2 L) y_{t+1} = u_t \\tag{6}\n",
    "$$\n",
    "\n",
    "where $ \\{u_t\\} $ is a bounded sequence, $ y_0 $ is an initial\n",
    "condition, $ | \\lambda_1 | < 1 $ and $ | \\lambda_2| >1 $.\n",
    "\n",
    "We seek a bounded sequence $ \\{y_t\\}_{t=0}^\\infty $ that satisfies\n",
    "[(6)](#equation-equn-6). Using insights from our analysis of the first-order equation,\n",
    "operate on both sides of [(6)](#equation-equn-6) by the forward inverse of\n",
    "$ (1-\\lambda_2 L) $ to rewrite equation [(6)](#equation-equn-6) as\n",
    "\n",
    "$$\n",
    "(1-\\lambda_1 L) y_{t+1} = -{\\frac{\\lambda_2^{-1}}{1 - \\lambda_2^{-1}L^{-1}}} u_{t+1}\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "\n",
    "<a id='equation-equn-7'></a>\n",
    "$$\n",
    "y_{t+1} = \\lambda_1 y_t - \\lambda_2^{-1} \\sum_{j=0}^\\infty \\lambda_2^{-j} u_{t+j+1} . \\tag{7}\n",
    "$$\n",
    "\n",
    "Thus, we obtained equation [(7)](#equation-equn-7) by\n",
    "solving stable roots (in this case $ \\lambda_1 $) **backward**, and\n",
    "unstable roots (in this case $ \\lambda_2 $) **forward**.\n",
    "\n",
    "Equation [(7)](#equation-equn-7) has a form that we shall encounter often.\n",
    "\n",
    "$ \\lambda_1 y_t $ is called the **feedback part** and\n",
    "$ -{\\frac{\\lambda_2^{-1}}{1 - \\lambda_2^{-1}L^{-1}}} u_{t+1} $ is\n",
    "called the **feedforward part** of the solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Illustration: Cagan’s Model\n",
    "\n",
    "Let\n",
    "\n",
    "- $ m_t^d $ be the log of the demand for money  \n",
    "- $ m_t $ be the log of the supply of money  \n",
    "- $ p_t $ be the log of the price level  \n",
    "\n",
    "\n",
    "It follows that $ p_{t+1} - p_t $ is the rate of inflation.\n",
    "\n",
    "The logarithm of the demand for real money balances $ m_t^d - p_t $\n",
    "is an inverse function of the expected rate of inflation\n",
    "$ p_{t+1} - p_t $ for $ t \\geq 0 $:\n",
    "\n",
    "$$\n",
    "m_t^d - p_t = - \\beta (p_{t+1} - p_t ), \\quad \\beta >0\n",
    "$$\n",
    "\n",
    "Equate the demand for log money $ m_t^d $ to the supply of log money\n",
    "$ m_t $ in the above equation and rearrange to deduce that the\n",
    "logarithm of the price level $ p_t $ is related to the logarithm of\n",
    "the money supply $ m_t $ by\n",
    "\n",
    "\n",
    "<a id='equation-equation-1'></a>\n",
    "$$\n",
    "p_t = (1 -\\lambda) m_t + \\lambda p_{t+1} \\tag{8}\n",
    "$$\n",
    "\n",
    "where $ \\lambda \\equiv \\frac{\\beta}{1+\\beta} \\in (0,1) $.\n",
    "\n",
    "Solving the first order difference equation [(8)](#equation-equation-1) forward gives\n",
    "\n",
    "\n",
    "<a id='equation-equation-2'></a>\n",
    "$$\n",
    "p_t = (1 - \\lambda) \\sum_{j=0}^\\infty \\lambda^j m_{t+j}, \\tag{9}\n",
    "$$\n",
    "\n",
    "which is the unique **stable** solution of difference equation [(8)](#equation-equation-1) among\n",
    "a class of more general solutions\n",
    "\n",
    "\n",
    "<a id='equation-equation-1a'></a>\n",
    "$$\n",
    "p_t = (1 - \\lambda) \\sum_{j=0}^\\infty \\lambda^j m_{t+j} + c \\lambda^{-t} \\tag{10}\n",
    "$$\n",
    "\n",
    "that is indexed by the real number $ c \\in {\\bf R} $.\n",
    "\n",
    "Because we want to focus on stable solutions, we set $ c=0 $.\n",
    "\n",
    "We begin by assuming that the log of the money supply is **exogenous**\n",
    "in the sense that it is an autonomous process that does not feed back on\n",
    "the log of the price level.\n",
    "\n",
    "In particular, we assume that the log of the money supply is described\n",
    "by the linear state space system\n",
    "\n",
    "\n",
    "<a id='equation-equation-3'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "  m_t &  = G x_t \\\\ x_{t+1} & = A x_t\n",
    " \\end{aligned} \\tag{11}\n",
    "$$\n",
    "\n",
    "where $ x_t $ is an $ n \\times 1 $ vector that does not include\n",
    "$ p_t $ or lags of $ p_t $, $ A $ is an $ n \\times n $\n",
    "matrix with eigenvalues that are less than $ \\lambda^{-1} $ in\n",
    "absolute values, and $ G $ is a $ 1 \\times n $ selector matrix.\n",
    "\n",
    "Variables appearing in the vector $ x_t $ contain information that\n",
    "might help predict future values of the money supply.\n",
    "\n",
    "We’ll take an example in which $ x_t $ includes only $ m_t $,\n",
    "possibly lagged values of $ m $, and a constant.\n",
    "\n",
    "An example of such an $ \\{m_t\\} $ process that fits info state space\n",
    "system [(11)](#equation-equation-3) is one that satisfies the second order linear difference\n",
    "equation\n",
    "\n",
    "$$\n",
    "m_{t+1} = \\alpha + \\rho_1 m_t + \\rho_2 m_{t-1}\n",
    "$$\n",
    "\n",
    "where the zeros of the characteristic polynomial\n",
    "$ (1 - \\rho_1 z - \\rho_2 z^2) $ are strictly greater than $ 1 $\n",
    "in modulus\n",
    "\n",
    "We seek a stable or non-explosive solution of the difference equation [(8)](#equation-equation-1) that\n",
    "obeys the system comprised of [(8)](#equation-equation-1)-[(11)](#equation-equation-3).\n",
    "\n",
    "By stable or non-explosive, we mean that neither $ m_t $ nor $ p_t $\n",
    "diverges as $ t \\rightarrow + \\infty $.\n",
    "\n",
    "This means that we are shutting down the term $ c \\lambda^{-t} $ in equation [(10)](#equation-equation-1a) above by setting $ c=0 $\n",
    "\n",
    "The solution we are after is\n",
    "\n",
    "\n",
    "<a id='equation-equation-4'></a>\n",
    "$$\n",
    "p_t = F x_t \\tag{12}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "\n",
    "<a id='equation-equation-5'></a>\n",
    "$$\n",
    "F = (1-\\lambda) G (I - \\lambda A)^{-1} \\tag{13}\n",
    "$$\n",
    "\n",
    "**Note:** As mentioned above, an *explosive solution* of difference\n",
    "equation [(8)](#equation-equation-1) can be constructed by adding to the right hand of [(12)](#equation-equation-4) a\n",
    "sequence $ c \\lambda^{-t} $ where $ c $ is an arbitrary positive\n",
    "constant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some Python code\n",
    "\n",
    "We’ll construct examples that illustrate [(11)](#equation-equation-3).\n",
    "\n",
    "Our first example takes as the law of motion for the log money supply\n",
    "the second order difference equation\n",
    "\n",
    "\n",
    "<a id='equation-equation-6'></a>\n",
    "$$\n",
    "m_{t+1} = \\alpha + \\rho_1 m_t + \\rho_2 m_{t-1} \\tag{14}\n",
    "$$\n",
    "\n",
    "that is parameterized by $ \\rho_1, \\rho_2, \\alpha $\n",
    "\n",
    "To capture this parameterization with system [(9)](#equation-equation-2) we set\n",
    "\n",
    "$$\n",
    "x_t = \\begin{bmatrix} 1 \\cr m_t \\cr m_{t-1} \\end{bmatrix} , \\quad\n",
    "  A= \\begin{bmatrix} 1 & 0 & 0 \\cr\n",
    "                     \\alpha & \\rho_1 & \\rho_2 \\cr\n",
    "                      0 & 1 & 0 \\end{bmatrix} , \\quad\n",
    "  G = \\begin{bmatrix} 0 & 1 & 0 \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Here is Python code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "λ = .9\n",
    "\n",
    "α = 0\n",
    "ρ1 = .9\n",
    "ρ2 = .05\n",
    "\n",
    "A = np.array([[1,  0,  0],\n",
    "              [α, ρ1, ρ2],\n",
    "              [0,  1,  0]])\n",
    "G = np.array([[0, 1, 0]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The matrix $ A $ has one eigenvalue equal to unity that is\n",
    "associated with the $ A_{11} $ component that captures a\n",
    "constant component of the state $ x_t $.\n",
    "\n",
    "We can verify that the two eigenvalues of $ A $ not associated with\n",
    "the constant in the state $ x_t $ are strictly less than unity in\n",
    "modulus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "eigvals = np.linalg.eigvals(A)\n",
    "print(eigvals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "(abs(eigvals) <= 1).all()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let’s compute $ F $ in formulas [(12)](#equation-equation-4) and [(13)](#equation-equation-5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# compute the solution, i.e. forumula (3)\n",
    "F = (1 - λ) * G @ np.linalg.inv(np.eye(A.shape[0]) - λ * A)\n",
    "print(\"F= \",F)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let’s simulate paths of $ m_t $ and $ p_t $ starting from an\n",
    "initial value $ x_0 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# set the initial state\n",
    "x0 = np.array([1, 1, 0])\n",
    "\n",
    "T = 100 # length of simulation\n",
    "\n",
    "m_seq = np.empty(T+1)\n",
    "p_seq = np.empty(T+1)\n",
    "\n",
    "m_seq[0] = G @ x0\n",
    "p_seq[0] = F @ x0\n",
    "\n",
    "# simulate for T periods\n",
    "x_old = x0\n",
    "for t in range(T):\n",
    "\n",
    "    x = A @ x_old\n",
    "\n",
    "    m_seq[t+1] = G @ x\n",
    "    p_seq[t+1] = F @ x\n",
    "\n",
    "    x_old = x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(range(T+1), m_seq, label='$m_t$')\n",
    "plt.plot(range(T+1), p_seq, label='$p_t$')\n",
    "plt.xlabel('t')\n",
    "plt.title(f'λ={λ}, α={α}, $ρ_1$={ρ1}, $ρ_2$={ρ2}')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above graph, why is the log of the price level always less than\n",
    "the log of the money supply?\n",
    "\n",
    "The answer is because\n",
    "\n",
    "- according to equation [(9)](#equation-equation-2), $ p_t $ is a geometric weighted\n",
    "  average of current and future values of $ m_t $, and  \n",
    "- it happens that in this example future $ m $’s are always less\n",
    "  than the current $ m $  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alternative code\n",
    "\n",
    "We could also have run the simulation using the quantecon\n",
    "**LinearStateSpace** code.\n",
    "\n",
    "The following code block performs the calculation with that code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# construct a LinearStateSpace instance\n",
    "\n",
    "# stack G and F\n",
    "G_ext = np.vstack([G, F])\n",
    "\n",
    "C = np.zeros((A.shape[0], 1))\n",
    "\n",
    "ss = qe.LinearStateSpace(A, C, G_ext, mu_0=x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "T = 100\n",
    "\n",
    "# simulate using LinearStateSpace\n",
    "x, y = ss.simulate(ts_length=T)\n",
    "\n",
    "# plot\n",
    "plt.figure()\n",
    "plt.plot(range(T), y[0,:], label='$m_t$')\n",
    "plt.plot(range(T), y[1,:], label='$p_t$')\n",
    "plt.xlabel('t')\n",
    "plt.title(f'λ={λ}, α={α}, $ρ_1$={ρ1}, $ρ_2$={ρ2}')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Special case\n",
    "\n",
    "To simplify our presentation in ways that will let focus on an important\n",
    "idea, in the above second-order difference equation [(14)](#equation-equation-6) that governs\n",
    "$ m_t $, we now set $ \\alpha =0 $,\n",
    "$ \\rho_1 = \\rho \\in (-1,1) $, and $ \\rho_2 =0 $ so that the law\n",
    "of motion for $ m_t $ becomes\n",
    "\n",
    "\n",
    "<a id='equation-equation-7'></a>\n",
    "$$\n",
    "m_{t+1} =\\rho m_t \\tag{15}\n",
    "$$\n",
    "\n",
    "and the state $ x_t $ becomes\n",
    "\n",
    "$$\n",
    "x_t = m_t .\n",
    "$$\n",
    "\n",
    "Consequently,  we can set $ G =1, A =\\rho $ making our formula [(13)](#equation-equation-5) for $ F $\n",
    "become\n",
    "\n",
    "$$\n",
    "F = (1-\\lambda) (1 -\\lambda \\rho)^{-1} .\n",
    "$$\n",
    "\n",
    "and the log the log price level satisfies\n",
    "\n",
    "$$\n",
    "p_t = F m_t .\n",
    "$$\n",
    "\n",
    "Please keep these formulas in mind as we investigate an alternative\n",
    "route to and interpretation of the formula for $ F $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Another perspective\n",
    "\n",
    "Above, we imposed stability or non-explosiveness on the solution of the key difference equation [(8)](#equation-equation-1)\n",
    "in Cagan’s model by solving the  unstable root $ \\lambda^{-1} $ forward.\n",
    "\n",
    "To shed light on the mechanics involved in imposing stability on a\n",
    "solution of a potentially unstable system of linear difference equations\n",
    "and to prepare the way for generalizations of our model in which the\n",
    "money supply is allowed to feed back on the price level itself, we stack\n",
    "equations [(8)](#equation-equation-1) and [(15)](#equation-equation-7) to form the system\n",
    "\n",
    "\n",
    "<a id='equation-equation-8'></a>\n",
    "$$\n",
    "\\begin{bmatrix} m_{t+1} \\cr p_{t+1} \\end{bmatrix} = \\begin{bmatrix} \\rho & 0 \\\\ - (1-\\lambda)/\\lambda & \\lambda^{-1}  \\end{bmatrix} \\begin{bmatrix} m_t \\\\ p_t \\end{bmatrix} \\tag{16}\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "\n",
    "<a id='equation-equation-9'></a>\n",
    "$$\n",
    "y_{t+1} = H y_t, \\quad t \\geq 0 \\tag{17}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "\n",
    "<a id='equation-equation-10'></a>\n",
    "$$\n",
    "H = \\begin{bmatrix} \\rho & 0 \\\\ - (1-\\lambda)/\\lambda & \\lambda^{-1}  \\end{bmatrix} . \\tag{18}\n",
    "$$\n",
    "\n",
    "Transition matrix $ H $ has eigenvalues $ \\rho \\in (0,1) $ and\n",
    "$ \\lambda^{-1} > 1 $.\n",
    "\n",
    "Because an eigenvalue of $ H $ exceeds unity, if we iterate on\n",
    "equation [(17)](#equation-equation-9) starting from an arbitrary initial vector\n",
    "$ y_0 = \\begin{bmatrix} m_0 \\\\ p_0 \\end{bmatrix} $, we discover that\n",
    "in general absolute values of both components of $ y_t $ diverge\n",
    "toward $ +\\infty $ as $ t \\rightarrow + \\infty $.\n",
    "\n",
    "To substantiate this claim, we can use the eigenector matrix\n",
    "decomposition of $ H $ that is available to us because the\n",
    "eigenvalues of $ H $ are distinct\n",
    "\n",
    "$$\n",
    "H = Q \\Lambda Q^{-1} .\n",
    "$$\n",
    "\n",
    "Here $ \\Lambda $ is a diagonal matrix of eigenvalues of $ H $\n",
    "and $ Q $ is a matrix whose columns are eigenvectors of the\n",
    "corresponding eigenvalues.\n",
    "\n",
    "Note that\n",
    "\n",
    "$$\n",
    "H^t = Q \\Lambda^t Q^{-1}\n",
    "$$\n",
    "\n",
    "so that\n",
    "\n",
    "$$\n",
    "y_t = Q \\Lambda^t Q^{-1} y_0\n",
    "$$\n",
    "\n",
    "For almost all initial vectors $ y_0 $, the presence of the\n",
    "eigenvalue $ \\lambda^{-1} > 1 $ causes both components of\n",
    "$ y_t $ to diverge in absolute value to $ +\\infty $.\n",
    "\n",
    "To explore this outcome in more detail, we use the following\n",
    "transformation\n",
    "\n",
    "$$\n",
    "y^*_t = Q^{-1} y_t\n",
    "$$\n",
    "\n",
    "that allows us to represent the dynamics in a way that isolates the\n",
    "source of the propensity of paths to diverge:\n",
    "\n",
    "$$\n",
    "y^*_{t+1} = \\Lambda^t y^*_t\n",
    "$$\n",
    "\n",
    "Staring at this equation indicates that unless\n",
    "\n",
    "\n",
    "<a id='equation-equation-11'></a>\n",
    "$$\n",
    "y^*_0 = \\begin{bmatrix} y^*_{1,0} \\cr 0 \\end{bmatrix} , \\tag{19}\n",
    "$$\n",
    "\n",
    "the path of $ y^*_t $ and therefore the paths of both components of\n",
    "$ y_t = Q y^*_t $ will diverge in absolute value as\n",
    "$ t \\rightarrow +\\infty $. (We say that the paths *explode*)\n",
    "\n",
    "Equation [(19)](#equation-equation-11) also leads us to conclude that there is a unique setting\n",
    "for the initial vector $ y_0 $ for which both components of\n",
    "$ y_t $ do not diverge.\n",
    "\n",
    "The required setting of $ y_0 $ must evidently have the property\n",
    "that\n",
    "\n",
    "$$\n",
    "Q y_0 =  y^*_0 = \\begin{bmatrix} y^*_{1,0} \\cr 0 \\end{bmatrix} .\n",
    "$$\n",
    "\n",
    "But note that since\n",
    "$ y_0 = \\begin{bmatrix} m_0 \\cr p_0 \\end{bmatrix} $ and $ m_0 $\n",
    "is given to us an an initial condition, it has to be $ p_0 $ that\n",
    "does all the adjusting to satisfy this equation.\n",
    "\n",
    "Sometimes this situation is described by saying that while $ m_0 $\n",
    "is truly a **state** variable, $ p_0 $ is a **jump** variable that\n",
    "is free to adjust at $ t=0 $ in order to satisfy the equation.\n",
    "\n",
    "Thus, in a nutshell the unique value of the vector $ y_0 $ for which\n",
    "the paths of $ y_t $ do not diverge must have second component\n",
    "$ p_0 $ that verifies equality [(19)](#equation-equation-11) by setting the second component\n",
    "of $ y^*_0 $ equal to zero.\n",
    "\n",
    "The component $ p_0 $ of the initial vector\n",
    "$ y_0 = \\begin{bmatrix} m_0 \\cr p_0 \\end{bmatrix} $ must evidently\n",
    "satisfy\n",
    "\n",
    "$$\n",
    "Q^{\\{2\\}} y_0 =0\n",
    "$$\n",
    "\n",
    "where $ Q^{\\{2\\}} $ denotes the second row of $ Q^{-1} $, a\n",
    "restriction that is equivalent to\n",
    "\n",
    "\n",
    "<a id='equation-equation-12'></a>\n",
    "$$\n",
    "Q^{21} m_0 + Q^{22} p_0 = 0 \\tag{20}\n",
    "$$\n",
    "\n",
    "where $ Q^{ij} $ denotes the $ (i,j) $ component of\n",
    "$ Q^{-1} $.\n",
    "\n",
    "Solving this equation for $ p_0 $ we find\n",
    "\n",
    "\n",
    "<a id='equation-equation-13'></a>\n",
    "$$\n",
    "p_0 = - (Q^{22})^{-1} Q^{21} m_0. \\tag{21}\n",
    "$$\n",
    "\n",
    "This is the unique **stabilizing value** of $ p_0 $ as a function of\n",
    "$ m_0 $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Refining the formula\n",
    "\n",
    "We can get an even more convenient formula for $ p_0 $ that is cast\n",
    "in terms of components of $ Q $ instead of components of\n",
    "$ Q^{-1} $.\n",
    "\n",
    "To get this formula, first note that because $ (Q^{21}\\ Q^{22}) $ is\n",
    "the second row of the inverse of $ Q $ and because\n",
    "$ Q^{-1} Q = I $, it follows that\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix} Q^{21} & Q^{22} \\end{bmatrix}  \\begin{bmatrix} Q_{11}\\cr Q_{21} \\end{bmatrix} = 0\n",
    "$$\n",
    "\n",
    "which implies that\n",
    "\n",
    "$$\n",
    "Q^{21} Q_{11} + Q^{22} Q_{21} = 0.\n",
    "$$\n",
    "\n",
    "Therefore,\n",
    "\n",
    "$$\n",
    "-(Q^{22})^{-1} Q^{21} = Q_{21} Q^{-1}_{11}.\n",
    "$$\n",
    "\n",
    "So we can write\n",
    "\n",
    "\n",
    "<a id='equation-equation-14'></a>\n",
    "$$\n",
    "p_0 = Q_{21} Q_{11}^{-1} m_0 . \\tag{22}\n",
    "$$\n",
    "\n",
    "It can be verified that this formula replicates itself over time so that\n",
    "\n",
    "\n",
    "<a id='equation-equation-15'></a>\n",
    "$$\n",
    "p_t = Q_{21} Q^{-1}_{11} m_t. \\tag{23}\n",
    "$$\n",
    "\n",
    "To implement formula [(23)](#equation-equation-15), we want to compute $ Q_1 $ the\n",
    "eigenvector of $ Q $ associated with the stable eigenvalue\n",
    "$ \\rho $ of $ Q $.\n",
    "\n",
    "By hand it can be verified that the eigenvector associated with the\n",
    "stable eigenvalue $ \\rho $ is proportional to\n",
    "\n",
    "$$\n",
    "Q_1  = \\begin{bmatrix} 1-\\lambda  \\rho \\\\ 1 - \\lambda   \\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "Notice that if we set $ A=\\rho $ and $ G=1 $ in our earlier\n",
    "formula for $ p_t $ we get\n",
    "\n",
    "$$\n",
    "Q = G (I - \\lambda A)^{-1} m_t =  (1-\\lambda) (1 - \\lambda \\rho)^{-1} m_t\n",
    "$$\n",
    "\n",
    "a formula that is equivalent with\n",
    "\n",
    "$$\n",
    "p_t = Q_{21} Q_{11}^{-1}  m_t ,\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "Q_1 = \\begin{bmatrix} Q_{11} \\\\ Q_{21}  \\end{bmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Some remarks about feedback\n",
    "\n",
    "We have expressed [(16)](#equation-equation-8) in what superficially appears to be a form in\n",
    "which $ y_{t+1} $ feeds back on $ y_t $. even though what we\n",
    "actually want to represent is that the component $ p_t $ feeds\n",
    "**forward** on $ p_{t+1} $, and through it, on future\n",
    "$ m_{t+j} $, $ j = 0, 1, 2, \\ldots $.\n",
    "\n",
    "A tell-tale sign that we should look beyond its superficial “feedback”\n",
    "form is that $ \\lambda^{-1} > 1 $ so that the matrix $ H $ in\n",
    "[(16)](#equation-equation-8) is **unstable**\n",
    "\n",
    "- it has one eigenvalue $ \\rho $ that is less than one in modulus\n",
    "  that does not imperil stability, but $ \\ldots $  \n",
    "- it has a second eigenvalue $ \\lambda^{-1} $ that exceeds one in\n",
    "  modulus and that makes $ H $ an unstable matrix  \n",
    "\n",
    "\n",
    "We’ll keep these observations in mind as we turn now to a case in which\n",
    "the log money supply actually does feed back on the log of the price\n",
    "level."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Log money supply feeds back on log price level\n",
    "\n",
    "The same pattern of eigenvalues splitting around unity, with one being\n",
    "below unity and another greater than unity, sometimes continues to\n",
    "prevail when there is  *feedback* from the log price level to the log\n",
    "money supply.\n",
    "\n",
    "Let the feedback rule be\n",
    "\n",
    "\n",
    "<a id='equation-equation-16'></a>\n",
    "$$\n",
    "m_{t+1} =  \\rho m_t + \\delta p_t \\tag{24}\n",
    "$$\n",
    "\n",
    "where $ \\rho \\in (0,1) $ as before and where we shall now allow\n",
    "$ \\delta \\neq 0 $.\n",
    "\n",
    "However,\n",
    "$ \\delta $ cannot be too large if things are to fit together as we\n",
    "wish to deliver a stable system for some initial value $ p_0 $ that we want to determine uniquely.\n",
    ".\n",
    "\n",
    "The forward-looking equation [(8)](#equation-equation-1) continues to describe equality between\n",
    "the demand and supply of money.\n",
    "\n",
    "We assume that equations [(8)](#equation-equation-1) and [(24)](#equation-equation-16) govern\n",
    "$ y_t \\equiv \\begin{bmatrix} m_t \\cr p_t \\end{bmatrix} $ for\n",
    "$ t \\geq 0 $\n",
    "\n",
    "The transition matrix $ H $ in the law of motion\n",
    "\n",
    "$$\n",
    "y_{t+1} = H y_t\n",
    "$$\n",
    "\n",
    "now becomes\n",
    "\n",
    "$$\n",
    "H = \\begin{bmatrix} \\rho & \\delta \\\\ - (1-\\lambda)/\\lambda & \\lambda^{-1}  \\end{bmatrix} .\n",
    "$$\n",
    "\n",
    "We take $ m_0 $ as a given intial condition and as before seek an\n",
    "initial value $ p_0 $ that stabilizes the system in the sense that\n",
    "$ y_t $ converges as $ t \\rightarrow + \\infty $.\n",
    "\n",
    "Our approach is identical with that followed above and is based on an\n",
    "eigenvalue decomposition in which, cross our fingers, one eigenvalue\n",
    "exceeds unity and the other is less than unity in absolute value.\n",
    "\n",
    "When $ \\delta \\neq 0 $ as we now assume, the eigenvalues of\n",
    "$ H $ are no longer $ \\rho \\in (0,1) $ and\n",
    "$ \\lambda^{-1} > 1 $\n",
    "\n",
    "We’ll just calculate them and apply the same algorithm that we used\n",
    "above.\n",
    "\n",
    "That algorithm remains valid so long as the eigenvalues split around\n",
    "unity as before.\n",
    "\n",
    "Again we assume that $ m_0 $ is  an initial condition, but that\n",
    "$ p_0 $ is not given but to be solved for.\n",
    "\n",
    "Let’s write and execute some Python code that will let us explore how outcomes depend on\n",
    "$ \\delta $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def construct_H(ρ, λ, δ):\n",
    "    \"contruct matrix H given parameters.\"\n",
    "\n",
    "    H = np.empty((2, 2))\n",
    "    H[0, :] = ρ,δ\n",
    "    H[1, :] = - (1 - λ) / λ, 1 / λ\n",
    "\n",
    "    return H\n",
    "\n",
    "def H_eigvals(ρ=.9, λ=.5, δ=0):\n",
    "    \"compute the eigenvalues of matrix H given parameters.\"\n",
    "\n",
    "    # construct H matrix\n",
    "    H = construct_H(ρ, λ, δ)\n",
    "\n",
    "    # compute eigenvalues\n",
    "    eigvals = np.linalg.eigvals(H)\n",
    "\n",
    "    return eigvals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H_eigvals()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that a negative δ will not imperil the stability of the matrix\n",
    "$ H $, even if it has a big absolute value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# small negative δ\n",
    "H_eigvals(δ=-0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# large negative δ\n",
    "H_eigvals(δ=-1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A sufficiently small positive δ also causes no problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# sufficiently small positive δ\n",
    "H_eigvals(δ=0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But a large enough positive δ makes both eigenvalues of $ H $\n",
    "strictly greater than unity in modulus.\n",
    "\n",
    "For example,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H_eigvals(δ=0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to study systems in which one eigenvalue exceeds unity in\n",
    "modulus while the other is less than unity in modulus, so we avoid\n",
    "values of $ \\delta $ that are too large"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def magic_p0(m0, ρ=.9, λ=.5, δ=0):\n",
    "    \"\"\"\n",
    "    Use the magic formula (8) to compute the level of p0\n",
    "    that makes the system stable.\n",
    "    \"\"\"\n",
    "\n",
    "    H = construct_H(ρ, λ, δ)\n",
    "    eigvals, Q = np.linalg.eig(H)\n",
    "\n",
    "    # find the index of the smaller eigenvalue\n",
    "    ind = 0 if eigvals[0] < eigvals[1] else 1\n",
    "\n",
    "    # verify that the eigenvalue is less than unity\n",
    "    if eigvals[ind] > 1:\n",
    "\n",
    "        print(\"both eigenvalues exceed unity in modulus\")\n",
    "\n",
    "        return None\n",
    "\n",
    "    p0 = Q[1, ind] / Q[0, ind] * m0\n",
    "\n",
    "    return p0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s plot how the solution $ p_0 $ changes as $ m_0 $\n",
    "changes for different settings of $ \\delta $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "m_range = np.arange(0.1, 2., 0.1)\n",
    "\n",
    "for δ in [-0.05, 0, 0.05]:\n",
    "    plt.plot(m_range, [magic_p0(m0, δ=δ) for m0 in m_range], label=f\"δ={δ}\")\n",
    "plt.legend()\n",
    "\n",
    "plt.xlabel(\"$m_0$\")\n",
    "plt.ylabel(\"$p_0$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To look at things from a different angle, we can fix the initial value $ m_0 $ and\n",
    "see how $ p_0 $ changes as $ \\delta $ changes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "m0 = 1\n",
    "\n",
    "δ_range = np.linspace(-0.05, 0.05, 100)\n",
    "plt.plot(δ_range, [magic_p0(m0, δ=δ) for δ in δ_range])\n",
    "plt.xlabel('$\\delta$')\n",
    "plt.ylabel('$p_0$')\n",
    "plt.title(f'$m_0$={m0}')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that when $ \\delta $ is large enough, both eigenvalues exceed\n",
    "unity in modulus, causing a stabilizing value of $ p_0 $ not to\n",
    "exist."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "magic_p0(1, δ=0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Big $ P $, little $ p $ interpretation\n",
    "\n",
    "It is helpful to view our solutions with feedback from the price level or inflation to money or the rate of money\n",
    "creation in terms of the Big $ K $, little $ k $ idea discussed in [Rational Expectations Models](https://python-programming.quantecon.org/rational_expectations.html)\n",
    "\n",
    "This will help us sort out what is taken as given by the decision makers who use the\n",
    "difference equation [(9)](#equation-equation-2) to determine $ p_t $ as a function of their forecasts of future values of\n",
    "$ m_t $.\n",
    "\n",
    "Let’s write the stabilizing solution that we have computed using the eigenvector decomposition of $ H $ as\n",
    "$ P_t = F^* m_t $ where\n",
    "\n",
    "$$\n",
    "F^* = Q_{21} Q_{11}^{-1}\n",
    "$$\n",
    "\n",
    "Then from $ P_{t+1} = F^* m_{t+1} $ and $ m_{t+1} = \\rho m_t + \\delta P_t $ we can deduce the recursion $ P_{t+1} = F^* \\rho m_t + F^* \\delta P_t $ and create the stacked system\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix} m_{t+1} \\cr P_{t+1} \\end{bmatrix}  =    \\begin{bmatrix} \\rho & \\delta \\cr\n",
    "                F^* \\rho & F^* \\delta   \\end{bmatrix} \\begin{bmatrix} m_t \\cr P_t \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "$$\n",
    "x_{t+1} = A x_t\n",
    "$$\n",
    "\n",
    "where $ x_t = \\begin{bmatrix} m_t \\cr P_t \\end{bmatrix} $.\n",
    "\n",
    "Then apply formula [(13)](#equation-equation-5) for $ F $ to deduce that\n",
    "\n",
    "$$\n",
    "p_t = F \\begin{bmatrix} m_t \\cr P_t \\end{bmatrix} = F \\begin{bmatrix} m_t \\cr F^* m_t \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "which implies that\n",
    "\n",
    "$$\n",
    "p_t = \\begin{bmatrix} F_1 & F_2 \\end{bmatrix}    \\begin{bmatrix} m_t \\cr F^* m_t \\end{bmatrix} = F_1 m_t + F_2 F^* m_t\n",
    "$$\n",
    "\n",
    "so that we expect to have\n",
    "\n",
    "$$\n",
    "F^* = F_1 + F_2 F^*\n",
    "$$\n",
    "\n",
    "We verify this equality in the next block of Python code that implements the following\n",
    "computations.\n",
    "\n",
    "1. For the system with $ \\delta\\neq 0 $ so that there is feedback,\n",
    "  we compute the stabilizing solution for $ p_t $ in the form\n",
    "  $ p_t = F^* m_t $ where $ F^* = Q_{21}Q_{11}^{-1} $ as above.  \n",
    "1. Recalling the system [(11)](#equation-equation-3), [(12)](#equation-equation-4), and [(13)](#equation-equation-5) above, we define\n",
    "  $ x_t = \\begin{bmatrix} m_t \\cr P_t \\end{bmatrix} $ and notice\n",
    "  that it is Big $ P_t $ and not little $ p_t $ here. Then we form $ A $ and $ G $ as\n",
    "  $ A = \\begin{bmatrix}\\rho & \\delta \\cr F^* \\rho & F^*\\delta \\end{bmatrix} $\n",
    "  and $ G = \\begin{bmatrix} 1 & 0 \\end{bmatrix} $ and we compute\n",
    "  $ \\begin{bmatrix}  F_1 &  F_2 \\end{bmatrix} \\equiv F $\n",
    "  from equation [(13)](#equation-equation-5) above.  \n",
    "1. We compute $ F_1 +  F_2 F^* $ and compare it\n",
    "  with $ F^* $ and verify equality.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# set parameters\n",
    "ρ = .9\n",
    "λ = .5\n",
    "δ = .05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# solve for F_star\n",
    "H = construct_H(ρ, λ, δ)\n",
    "eigvals, Q = np.linalg.eig(H)\n",
    "\n",
    "ind = 0 if eigvals[0] < eigvals[1] else 1\n",
    "F_star = Q[1, ind] / Q[0, ind]\n",
    "F_star"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# solve for F_check\n",
    "A = np.empty((2, 2))\n",
    "A[0, :] = ρ, δ\n",
    "A[1, :] = F_star * A[0, :]\n",
    "\n",
    "G = np.array([1, 0])\n",
    "\n",
    "F_check= (1 - λ) * G @ np.linalg.inv(np.eye(2) - λ * A)\n",
    "F_check"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare $ F^* $ with $ F_1 + F_2 F^* $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "F_check[0] + F_check[1] * F_star, F_star"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fun with Sympy code\n",
    "\n",
    "This section is a small gift for readers who have made it this far.\n",
    "\n",
    "It puts Sympy to work on our model.\n",
    "\n",
    "Thus, we  use Sympy to compute some of the key objects comprising the eigenvector decomposition of $ H $.\n",
    "\n",
    "$ H $ with nonzero $ \\delta $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "λ, δ, ρ = symbols('λ, δ, ρ')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H1 = Matrix([[ρ,δ], [- (1 - λ) / λ, λ ** -1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H1.eigenvals()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H1.eigenvects()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ H $ with $ \\delta $ being zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H2 = Matrix([[ρ,0], [- (1 - λ) / λ, λ ** -1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H2.eigenvals()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H2.eigenvects()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we do induce sympy to do the following fun things for us analytically:\n",
    "\n",
    "1. We compute the matrix $ Q $ whose first column is\n",
    "  the eigenvector associated with $ \\rho $. and whose second column\n",
    "  is the eigenvector associated with $ \\lambda^{-1} $.  \n",
    "1. We use sympy to compute the inverse $ Q^{-1} $ of $ Q $\n",
    "  (both in symbols).  \n",
    "1. We use sympy to compute $ Q_{21} Q_{11}^{-1} $ (in symbols).  \n",
    "1. Where $ Q^{ij} $ denotes the $ (i,j) $ component of\n",
    "  $ Q^{-1} $, weighted use sympy to compute\n",
    "  $ - (Q^{22})^{-1} Q^{21} $ (again in symbols)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# construct Q\n",
    "vec = []\n",
    "for i, (eigval, _, eigvec) in enumerate(H2.eigenvects()):\n",
    "\n",
    "    vec.append(eigvec[0])\n",
    "\n",
    "    if eigval == ρ:\n",
    "        ind = i\n",
    "\n",
    "Q = vec[ind].col_insert(1, vec[1-ind])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "Q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ Q^{-1} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "Q_inv = Q ** (-1)\n",
    "Q_inv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ Q_{21}Q_{11}^{-1} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "Q[1, 0] / Q[0, 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ −(Q^{22})^{−1}Q^{21} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "- Q_inv[1, 0] / Q_inv[1, 1]"
   ]
  }
 ],
 "metadata": {
  "date": 1591697902.7693224,
  "filename": "re_with_feedback.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Stability in Linear Rational Expectations Models"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}