re_with_feedback.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "34276cb2",
   "metadata": {},
   "source": [
    "\n",
    "<a id='re-with-feedback'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "298d4861",
   "metadata": {},
   "source": [
    "# Stability in Linear Rational Expectations Models\n",
    "\n",
    "\n",
    "<a id='index-0'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "803daeaa",
   "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](#Fun-with-SymPy)  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "521f7da2",
   "metadata": {},
   "source": [
    "In addition to what’s in Anaconda, this lecture deploys the following libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9cacd9c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "!pip install quantecon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10b47bee",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams[\"figure.figsize\"] = (11, 5)  #set default figure size\n",
    "import numpy as np\n",
    "import quantecon as qe\n",
    "from sympy import init_printing, symbols, Matrix\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d613f62",
   "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 [[Cagan, 1956](https://python.quantecon.org/zreferences.html#id92)] linking\n",
    "the price level to the money supply.\n",
    "\n",
    "Cagan did not use a rational expectations version of his model, but Sargent [[Sargent, 1977](https://python.quantecon.org/zreferences.html#id93)] did.\n",
    "\n",
    "We study a rational expectations version of this model because it is intrinsically interesting and  because it\n",
    "has a mathematical structure that\n",
    "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",
    "The model determines  the price level or rate of inflation as a function of  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) [[Blanchard and Kahn, 1980](https://python.quantecon.org/zreferences.html#id251)] and Whiteman (1983) [[Whiteman, 1983](https://python.quantecon.org/zreferences.html#id252)] to solve a class of linear rational expectations models  \n",
    "- how to use **SymPy** to get analytical formulas for some key objects comprising a rational expectations equilibrium  \n",
    "\n",
    "\n",
    "Matrix decompositions employed here are described in more depth in this lecture [Lagrangian formulations](https://python.quantecon.org/lagrangian_lqdp.html).\n",
    "\n",
    "We formulate a version of  Cagan’s model under  rational expectations\n",
    "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",
   "id": "858180dc",
   "metadata": {},
   "source": [
    "## Linear Difference Equations\n",
    "\n",
    "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 deploy  $ L $  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",
    "We’ll often use the equality  $ L^{-1} x_t \\equiv x_{t+1} $ below.\n",
    "\n",
    "The algebra of lag and forward shift operators can simplify representing and solving linear difference equations."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5c738c2",
   "metadata": {},
   "source": [
    "### First Order\n",
    "\n",
    "We want to solve a linear first-order scalar difference equation.\n",
    "\n",
    "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{63.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{63.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 [(63.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 [(63.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{63.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{63.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{63.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 [(63.5)](#equation-equn-5) by $ (1 -\\lambda L) $ and compare to\n",
    "equation [(63.1)](#equation-equn-1).\n",
    "\n",
    "For any bounded $ \\{u_t\\} $ sequence, solution [(63.2)](#equation-equn-2) exists for $ |\\lambda | < 1 $ because\n",
    "the **distributed lag** in $ u $ converges.\n",
    "\n",
    "Solution [(63.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 [(63.2)](#equation-equn-2) may\n",
    "diverge, in which case a solution of this form does not exist.\n",
    "\n",
    "The distributed lead in $ u $ in [(63.5)](#equation-equn-5) need not\n",
    "converge when $ |\\lambda| < 1 $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d5af86c",
   "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{63.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",
    "[(63.6)](#equation-equn-6). Using insights from our analysis of the first-order equation,\n",
    "operate on both sides of [(63.6)](#equation-equn-6) by the forward inverse of\n",
    "$ (1-\\lambda_2 L) $ to rewrite equation [(63.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{63.7}\n",
    "$$\n",
    "\n",
    "Thus, we obtained equation [(63.7)](#equation-equn-7) by\n",
    "solving a stable root (in this case $ \\lambda_1 $) **backward**, and an\n",
    "unstable root (in this case $ \\lambda_2 $) **forward**.\n",
    "\n",
    "Equation [(63.7)](#equation-equn-7) has a form that we shall encounter often.\n",
    "\n",
    "- $ \\lambda_1 y_t $ is called the **feedback part**  \n",
    "- $ -{\\frac{\\lambda_2^{-1}}{1 - \\lambda_2^{-1}L^{-1}}} u_{t+1} $ is called the **feedforward part**  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98263a07",
   "metadata": {},
   "source": [
    "## Illustration: Cagan’s Model\n",
    "\n",
    "Now let’s use linear difference equations to represent and solve Sargent’s  [[Sargent, 1977](https://python.quantecon.org/zreferences.html#id93)] rational expectations version of\n",
    "Cagan’s model [[Cagan, 1956](https://python.quantecon.org/zreferences.html#id92)] that connects the price level to the public’s anticipations of future money supplies.\n",
    "\n",
    "Cagan did not use a rational expectations version of his model, but Sargent [[Sargent, 1977](https://python.quantecon.org/zreferences.html#id93)]\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{63.8}\n",
    "$$\n",
    "\n",
    "where $ \\lambda \\equiv \\frac{\\beta}{1+\\beta} \\in (0,1) $.\n",
    "\n",
    "(We note that the characteristic polynomial if $ 1 - \\lambda^{-1} z^{-1} = 0 $ so that the zero of the\n",
    "characteristic polynomial in this case is $ \\lambda \\in (0,1) $ which here is **inside** the unit circle.)\n",
    "\n",
    "Solving the first order difference equation [(63.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{63.9}\n",
    "$$\n",
    "\n",
    "which is the unique **stable** solution of difference equation [(63.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{63.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",
    "Equation [(63.10)](#equation-equation-1a) attributes  **perfect foresight** about the money supply sequence to the holders of real balances.\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{63.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 start with 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 [(63.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",
    "(Please see  [this](https://python.quantecon.org/samuelson.html) QuantEcon lecture for more about characteristic polynomials and their role in solving linear difference equations.)\n",
    "\n",
    "We seek a stable or non-explosive solution of the difference equation [(63.8)](#equation-equation-1) that\n",
    "obeys the system comprised of [(63.8)](#equation-equation-1)-[(63.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 requires  that we  shut down the term $ c \\lambda^{-t} $ in equation [(63.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{63.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{63.13}\n",
    "$$\n",
    "\n",
    ">**Note**\n",
    ">\n",
    ">As mentioned above, an *explosive solution* of difference\n",
    "equation [(63.8)](#equation-equation-1) can be constructed by adding to the right hand of [(63.12)](#equation-equation-4) a\n",
    "sequence $ c \\lambda^{-t} $ where $ c $ is an arbitrary positive\n",
    "constant."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37a1154b",
   "metadata": {},
   "source": [
    "## Some Python Code\n",
    "\n",
    "We’ll construct examples that illustrate [(63.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{63.14}\n",
    "$$\n",
    "\n",
    "that is parameterized by $ \\rho_1, \\rho_2, \\alpha $\n",
    "\n",
    "To capture this parameterization with system [(63.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,
   "id": "28fee5f1",
   "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",
   "id": "5ec2197b",
   "metadata": {},
   "source": [
    "The matrix $ A $ has one eigenvalue equal to unity.\n",
    "\n",
    "It is 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,
   "id": "1bc04070",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "eigvals = np.linalg.eigvals(A)\n",
    "print(eigvals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f84e801d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "(abs(eigvals) <= 1).all()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9628953d",
   "metadata": {},
   "source": [
    "Now let’s compute $ F $ in formulas [(63.12)](#equation-equation-4) and [(63.13)](#equation-equation-5)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "714e99cb",
   "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",
   "id": "7389c3a1",
   "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,
   "id": "630f33fd",
   "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,
   "id": "749d2a8c",
   "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",
   "id": "7d947725",
   "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",
    "Because\n",
    "\n",
    "- according to equation [(63.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",
   "id": "10b766f9",
   "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,
   "id": "619b3289",
   "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,
   "id": "e4431301",
   "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",
   "id": "f9da7957",
   "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 [(63.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{63.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 [(63.13)](#equation-equation-5) for $ F $\n",
    "become\n",
    "\n",
    "$$\n",
    "F = (1-\\lambda) (1 -\\lambda \\rho)^{-1} .\n",
    "$$\n",
    "\n",
    "so that 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 our formula for $ F $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34938f41",
   "metadata": {},
   "source": [
    "## Another Perspective\n",
    "\n",
    "Above, we imposed stability or non-explosiveness on the solution of the key difference equation [(63.8)](#equation-equation-1)\n",
    "in Cagan’s model by solving the  unstable root of the characteristic polynomial 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 [(63.8)](#equation-equation-1) and [(63.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{63.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{63.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{63.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 [(63.17)](#equation-equation-9) starting from an arbitrary initial vector\n",
    "$ y_0 = \\begin{bmatrix} m_0 \\\\ p_0 \\end{bmatrix} $ with $ m_0 >0, p_0 >0 $, 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 eigenvector 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 associated with 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 can 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{63.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 [(63.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 initial condition,  $ p_0 $ has to do 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",
    "must 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 [(63.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{63.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{63.21}\n",
    "$$\n",
    "\n",
    "This is the unique **stabilizing value** of $ p_0 $ expressed as a function of\n",
    "$ m_0 $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e43fb62c",
   "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{63.22}\n",
    "$$\n",
    "\n",
    "It can be verified that this formula replicates itself over time in the sense  that\n",
    "\n",
    "\n",
    "<a id='equation-equation-15'></a>\n",
    "$$\n",
    "p_t = Q_{21} Q^{-1}_{11} m_t. \\tag{63.23}\n",
    "$$\n",
    "\n",
    "To implement formula [(63.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",
    "p_t = 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",
   "id": "92a5aff4",
   "metadata": {},
   "source": [
    "### Remarks about Feedback\n",
    "\n",
    "We have expressed [(63.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",
    "[(63.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",
   "id": "3baffaf5",
   "metadata": {},
   "source": [
    "## Log money Supply Feeds Back on Log Price Level\n",
    "\n",
    "An arrangement of eigenvalues that split around unity, with one being\n",
    "below unity and another being greater than unity, sometimes prevails 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{63.24}\n",
    "$$\n",
    "\n",
    "where $ \\rho \\in (0,1) $  and where we shall now allow\n",
    "$ \\delta \\neq 0 $.\n",
    "\n",
    "**Warning:**  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, $ \\delta $ cannot be too large.\n",
    "\n",
    "The forward-looking equation [(63.8)](#equation-equation-1) continues to describe equality between\n",
    "the demand and supply of money.\n",
    "\n",
    "We assume that equations [(63.8)](#equation-equation-1) and [(63.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 initial 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 the one  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 $ will no longer be $ \\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,
   "id": "f67edc00",
   "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,
   "id": "1c5b200c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H_eigvals()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9b2495c",
   "metadata": {},
   "source": [
    "Notice that a negative $ \\delta $ will not imperil the stability of the matrix\n",
    "$ H $, even if it has a big absolute value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b6dc80f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# small negative δ\n",
    "H_eigvals(δ=-0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bab637e",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# large negative δ\n",
    "H_eigvals(δ=-1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "490ed1a4",
   "metadata": {},
   "source": [
    "A sufficiently small positive $ \\delta $ also causes no problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fb5a294",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# sufficiently small positive δ\n",
    "H_eigvals(δ=0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7963df15",
   "metadata": {},
   "source": [
    "But a large enough positive $ \\delta $ makes both eigenvalues of $ H $\n",
    "strictly greater than unity in modulus.\n",
    "\n",
    "For example,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2d1914b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H_eigvals(δ=0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "529a7a74",
   "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.\n",
    "\n",
    "That is, we want to avoid too much positive feedback from $ p_t $ to $ m_{t+1} $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b0965aa",
   "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",
   "id": "abcd92aa",
   "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,
   "id": "b35da828",
   "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",
   "id": "451be2d9",
   "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,
   "id": "8e6618a3",
   "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",
   "id": "a9860e07",
   "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,
   "id": "68175626",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "magic_p0(1, δ=0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "484e5c35",
   "metadata": {},
   "source": [
    "## Big $ P $, Little $ p $ Interpretation\n",
    "\n",
    "It is helpful to view our solutions of difference equations having  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.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 [(63.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",
    "Apply formula [(63.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 can anticipate that\n",
    "\n",
    "$$\n",
    "F^* = F_1 + F_2 F^*\n",
    "$$\n",
    "\n",
    "We shall 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 [(63.11)](#equation-equation-3), [(63.12)](#equation-equation-4), and [(63.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 [(63.13)](#equation-equation-5) above.  \n",
    "1. We compute $ F_1 +  F_2 F^* $ and compare it\n",
    "  with $ F^* $ and check for the anticipated  equality.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10b1d03e",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# set parameters\n",
    "ρ = .9\n",
    "λ = .5\n",
    "δ = .05"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "090b5824",
   "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,
   "id": "eadcd629",
   "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",
   "id": "da0abb32",
   "metadata": {},
   "source": [
    "Compare $ F^* $ with $ F_1 + F_2 F^* $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c94667b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "F_check[0] + F_check[1] * F_star, F_star"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eec38347",
   "metadata": {},
   "source": [
    "## Fun with SymPy\n",
    "\n",
    "This section is a  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  key objects comprising the eigenvector decomposition of $ H $.\n",
    "\n",
    "We start by generating an $ H $ with nonzero $ \\delta $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ea99e31",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "λ, δ, ρ = symbols('λ, δ, ρ')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "665ae4d6",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H1 = Matrix([[ρ,δ], [- (1 - λ) / λ, λ ** -1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a7a8552",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e062f01c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H1.eigenvals()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b15db31",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H1.eigenvects()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "400098b0",
   "metadata": {},
   "source": [
    "Now let’s compute $ H $ when  $ \\delta $ is zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bbc8ac5",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H2 = Matrix([[ρ,0], [- (1 - λ) / λ, λ ** -1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ea01d44",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccfb5125",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H2.eigenvals()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fe0f179",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "H2.eigenvects()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a16b5eee",
   "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} $, we use SymPy to compute\n",
    "  $ - (Q^{22})^{-1} Q^{21} $ (again in symbols)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "695710da",
   "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,
   "id": "ce1865d3",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "Q"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4003aeaa",
   "metadata": {},
   "source": [
    "$ Q^{-1} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6df6dad",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "Q_inv = Q ** (-1)\n",
    "Q_inv"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c8f2f6f",
   "metadata": {},
   "source": [
    "$ Q_{21}Q_{11}^{-1} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70f2662d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "Q[1, 0] / Q[0, 0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "717ca996",
   "metadata": {},
   "source": [
    "$ −(Q^{22})^{−1}Q^{21} $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2244a443",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "- Q_inv[1, 0] / Q_inv[1, 1]"
   ]
  }
 ],
 "metadata": {
  "date": 1714442509.9639528,
  "filename": "re_with_feedback.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Stability in Linear Rational Expectations Models"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}