lqcontrol.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "02b1f077",
   "metadata": {},
   "source": [
    "\n",
    "<a id='lqc'></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",
   "id": "bed629b0",
   "metadata": {},
   "source": [
    "# LQ Control: Foundations\n",
    "\n",
    "\n",
    "<a id='index-0'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01e84d1b",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [LQ Control: Foundations](#LQ-Control:-Foundations)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [Introduction](#Introduction)  \n",
    "  - [Optimality – Finite Horizon](#Optimality-–-Finite-Horizon)  \n",
    "  - [Implementation](#Implementation)  \n",
    "  - [Extensions and Comments](#Extensions-and-Comments)  \n",
    "  - [Further Applications](#Further-Applications)  \n",
    "  - [Exercises](#Exercises)  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bf680ec",
   "metadata": {},
   "source": [
    "In addition to what’s in Anaconda, this lecture will need the following libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f88009a",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "!pip install quantecon"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05087187",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "Linear quadratic (LQ) control refers to a class of dynamic optimization problems that have found applications in almost every scientific field.\n",
    "\n",
    "This lecture provides an introduction to LQ control and its economic applications.\n",
    "\n",
    "As we will see, LQ systems have a simple structure that makes them an excellent workhorse for a wide variety of economic problems.\n",
    "\n",
    "Moreover, while the linear-quadratic structure is restrictive, it is in fact far more flexible than it may appear initially.\n",
    "\n",
    "These themes appear repeatedly below.\n",
    "\n",
    "Mathematically, LQ control problems are closely related to [the Kalman filter](https://python.quantecon.org/kalman.html)\n",
    "\n",
    "- Recursive formulations of linear-quadratic control problems and Kalman filtering problems both involve matrix **Riccati equations**.  \n",
    "- Classical formulations of linear control and linear filtering problems make use of similar matrix decompositions (see for example [this lecture](https://python-advanced.quantecon.org/lu_tricks.html) and [this lecture](https://python-advanced.quantecon.org/classical_filtering.html)).  \n",
    "\n",
    "\n",
    "In reading what follows, it will be useful to have some familiarity with\n",
    "\n",
    "- matrix manipulations  \n",
    "- vectors of random variables  \n",
    "- dynamic programming and the Bellman equation (see for example [this lecture](https://intro.quantecon.org/short_path.html) and [this lecture](https://python.quantecon.org/optgrowth.html))  \n",
    "\n",
    "\n",
    "For additional reading on LQ control, see, for example,\n",
    "\n",
    "- [[Ljungqvist and Sargent, 2018](https://python.quantecon.org/zreferences.html#id186)], chapter 5  \n",
    "- [[Hansen and Sargent, 2008](https://python.quantecon.org/zreferences.html#id170)], chapter 4  \n",
    "- [[Hernandez-Lerma and Lasserre, 1996](https://python.quantecon.org/zreferences.html#id172)], section 3.5  \n",
    "\n",
    "\n",
    "In order to focus on computation, we leave longer proofs to these sources (while trying to provide as much intuition as possible).\n",
    "\n",
    "Let’s start with some imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c05153e",
   "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",
    "from quantecon import LQ"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3eb3255b",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "The “linear” part of LQ is a linear law of motion for the state, while the “quadratic” part refers to preferences.\n",
    "\n",
    "Let’s begin with the former, move on to the latter, and then put them together into an optimization problem."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dccd7a5",
   "metadata": {},
   "source": [
    "### The Law of Motion\n",
    "\n",
    "Let $ x_t $ be a vector describing the state of some economic system.\n",
    "\n",
    "Suppose that $ x_t $ follows a linear law of motion given by\n",
    "\n",
    "\n",
    "<a id='equation-lq-lom'></a>\n",
    "$$\n",
    "x_{t+1} = A x_t + B u_t + C w_{t+1},\n",
    "\\qquad t = 0, 1, 2, \\ldots \\tag{56.1}\n",
    "$$\n",
    "\n",
    "Here\n",
    "\n",
    "- $ u_t $ is a “control” vector, incorporating choices available to a decision-maker confronting the current state $ x_t $  \n",
    "- $ \\{w_t\\} $ is an uncorrelated zero mean shock process satisfying $ \\mathbb E w_t w_t' = I $, where the right-hand side is the identity matrix  \n",
    "\n",
    "\n",
    "Regarding the dimensions\n",
    "\n",
    "- $ x_t $ is $ n \\times 1 $, $ A $ is $ n \\times n $  \n",
    "- $ u_t $ is $ k \\times 1 $, $ B $ is $ n \\times k $  \n",
    "- $ w_t $ is $ j \\times 1 $, $ C $ is $ n \\times j $  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bd64d4c",
   "metadata": {},
   "source": [
    "#### Example 1\n",
    "\n",
    "Consider a household budget constraint given by\n",
    "\n",
    "$$\n",
    "a_{t+1} + c_t = (1 + r) a_t + y_t\n",
    "$$\n",
    "\n",
    "Here $ a_t $ is assets, $ r $ is a fixed interest rate, $ c_t $ is\n",
    "current consumption, and $ y_t $ is current non-financial income.\n",
    "\n",
    "If we suppose that $ \\{ y_t \\} $ is serially uncorrelated and $ N(0,\n",
    "\\sigma^2) $, then, taking $ \\{ w_t \\} $ to be standard normal, we can write\n",
    "the system as\n",
    "\n",
    "$$\n",
    "a_{t+1} = (1 + r) a_t - c_t + \\sigma w_{t+1}\n",
    "$$\n",
    "\n",
    "This is clearly a special case of [(56.1)](#equation-lq-lom), with assets being the state and consumption being the control.\n",
    "\n",
    "\n",
    "<a id='lq-hhp'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32f63ef4",
   "metadata": {},
   "source": [
    "#### Example 2\n",
    "\n",
    "One unrealistic feature of the previous model is that non-financial income has a zero mean and is often negative.\n",
    "\n",
    "This can easily be overcome by adding a sufficiently large mean.\n",
    "\n",
    "Hence in this example, we take $ y_t = \\sigma w_{t+1} + \\mu $ for some positive real number $ \\mu $.\n",
    "\n",
    "Another alteration that’s useful to introduce (we’ll see why soon) is to\n",
    "change the control variable from consumption\n",
    "to the deviation of consumption from some “ideal” quantity $ \\bar c $.\n",
    "\n",
    "(Most parameterizations will be such that $ \\bar c $ is large relative to the amount of consumption that is attainable in each period, and hence the household wants to increase consumption.)\n",
    "\n",
    "For this reason, we now take our control to be $ u_t := c_t - \\bar c $.\n",
    "\n",
    "In terms of these variables, the budget constraint $ a_{t+1} = (1 + r) a_t - c_t + y_t $ becomes\n",
    "\n",
    "\n",
    "<a id='equation-lq-lomwc'></a>\n",
    "$$\n",
    "a_{t+1} = (1 + r) a_t - u_t - \\bar c + \\sigma w_{t+1} + \\mu \\tag{56.2}\n",
    "$$\n",
    "\n",
    "How can we write this new system in the form of equation [(56.1)](#equation-lq-lom)?\n",
    "\n",
    "If, as in the previous example, we take $ a_t $ as the state, then we run into a problem:\n",
    "the law of motion contains some constant terms on the right-hand side.\n",
    "\n",
    "This means that we are dealing with an *affine* function, not a linear one\n",
    "(recall [this discussion](https://python.quantecon.org/linear_algebra.html#la-linear-map)).\n",
    "\n",
    "Fortunately, we can easily circumvent this problem by adding an extra state variable.\n",
    "\n",
    "In particular, if we write\n",
    "\n",
    "\n",
    "<a id='equation-lq-lowmc'></a>\n",
    "$$\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "a_{t+1} \\\\\n",
    "1\n",
    "\\end{array}\n",
    "\\right) =\n",
    "\\left(\n",
    "\\begin{array}{cc}\n",
    "1 + r & -\\bar c + \\mu \\\\\n",
    "0     & 1\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "a_t \\\\\n",
    "1\n",
    "\\end{array}\n",
    "\\right) +\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "-1 \\\\\n",
    "0\n",
    "\\end{array}\n",
    "\\right)\n",
    "u_t +\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "\\sigma \\\\\n",
    "0\n",
    "\\end{array}\n",
    "\\right)\n",
    "w_{t+1} \\tag{56.3}\n",
    "$$\n",
    "\n",
    "then the first row is equivalent to [(56.2)](#equation-lq-lomwc).\n",
    "\n",
    "Moreover, the model is now linear and can be written in the form of\n",
    "[(56.1)](#equation-lq-lom) by setting\n",
    "\n",
    "\n",
    "<a id='equation-lq-lowmc2'></a>\n",
    "$$\n",
    "x_t :=\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "a_t \\\\\n",
    "1\n",
    "\\end{array}\n",
    "\\right),\n",
    "\\quad\n",
    "A :=\n",
    "\\left(\n",
    "\\begin{array}{cc}\n",
    "1 + r & -\\bar c + \\mu \\\\\n",
    "0     & 1\n",
    "\\end{array}\n",
    "\\right),\n",
    "\\quad\n",
    "B :=\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "-1 \\\\\n",
    "0\n",
    "\\end{array}\n",
    "\\right),\n",
    "\\quad\n",
    "C :=\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "\\sigma \\\\\n",
    "0\n",
    "\\end{array}\n",
    "\\right) \\tag{56.4}\n",
    "$$\n",
    "\n",
    "In effect, we’ve bought ourselves linearity by adding another state."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae83ee43",
   "metadata": {},
   "source": [
    "### Preferences\n",
    "\n",
    "In the LQ model, the aim is to minimize flow of losses, where time-$ t $ loss is given by the quadratic expression\n",
    "\n",
    "\n",
    "<a id='equation-lq-pref-flow'></a>\n",
    "$$\n",
    "x_t' R x_t + u_t' Q u_t \\tag{56.5}\n",
    "$$\n",
    "\n",
    "Here\n",
    "\n",
    "- $ R $ is assumed to be $ n \\times n $, symmetric and nonnegative definite.  \n",
    "- $ Q $ is assumed to be $ k \\times k $, symmetric and positive definite.  \n",
    "\n",
    "\n",
    ">**Note**\n",
    ">\n",
    ">In fact, for many economic problems, the definiteness conditions on $ R $ and $ Q $ can be relaxed.  It is sufficient that certain submatrices of $ R $ and $ Q $ be nonnegative definite. See [[Hansen and Sargent, 2008](https://python.quantecon.org/zreferences.html#id170)] for details."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7aa25a8",
   "metadata": {},
   "source": [
    "#### Example 1\n",
    "\n",
    "A very simple example that satisfies these assumptions is to take $ R $\n",
    "and $ Q $ to be identity matrices so that current loss is\n",
    "\n",
    "$$\n",
    "x_t' I x_t + u_t' I u_t = \\| x_t \\|^2 + \\| u_t \\|^2\n",
    "$$\n",
    "\n",
    "Thus, for both the state and the control, loss is measured as squared distance from the origin.\n",
    "\n",
    "(In fact, the general case [(56.5)](#equation-lq-pref-flow) can also be understood in this\n",
    "way, but with $ R $ and $ Q $ identifying other – non-Euclidean – notions of “distance” from the zero vector.)\n",
    "\n",
    "Intuitively, we can often think of the state $ x_t $ as representing deviation from a target, such\n",
    "as\n",
    "\n",
    "- deviation of inflation from some target level  \n",
    "- deviation of a firm’s capital stock from some desired quantity  \n",
    "\n",
    "\n",
    "The aim is to put the state close to the target, while using  controls parsimoniously."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91001186",
   "metadata": {},
   "source": [
    "#### Example 2\n",
    "\n",
    "In the household problem [studied above](#lq-hhp), setting $ R=0 $\n",
    "and $ Q=1 $ yields preferences\n",
    "\n",
    "$$\n",
    "x_t' R x_t + u_t' Q u_t = u_t^2 = (c_t - \\bar c)^2\n",
    "$$\n",
    "\n",
    "Under this specification, the household’s current loss is the squared deviation of consumption from the ideal level $ \\bar c $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a40db2a",
   "metadata": {},
   "source": [
    "## Optimality – Finite Horizon\n",
    "\n",
    "\n",
    "<a id='index-1'></a>\n",
    "Let’s now be precise about the optimization problem we wish to consider, and look at how to solve it."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9493e244",
   "metadata": {},
   "source": [
    "### The Objective\n",
    "\n",
    "We will begin with the finite horizon case, with terminal time $ T \\in \\mathbb N $.\n",
    "\n",
    "In this case, the aim is to choose a sequence of controls $ \\{u_0, \\ldots, u_{T-1}\\} $ to minimize the objective\n",
    "\n",
    "\n",
    "<a id='equation-lq-object'></a>\n",
    "$$\n",
    "\\mathbb E \\,\n",
    "\\left\\{\n",
    "    \\sum_{t=0}^{T-1} \\beta^t (x_t' R x_t + u_t' Q u_t) + \\beta^T x_T' R_f x_T\n",
    "\\right\\} \\tag{56.6}\n",
    "$$\n",
    "\n",
    "subject to the law of motion [(56.1)](#equation-lq-lom) and initial state $ x_0 $.\n",
    "\n",
    "The new objects introduced here are $ \\beta $ and the matrix $ R_f $.\n",
    "\n",
    "The scalar $ \\beta $ is the discount factor, while $ x' R_f x $ gives terminal loss associated with state $ x $.\n",
    "\n",
    "Comments:\n",
    "\n",
    "- We assume $ R_f $ to be $ n \\times n $, symmetric and nonnegative definite.  \n",
    "- We allow $ \\beta = 1 $, and hence include the undiscounted case.  \n",
    "- $ x_0 $ may itself be random, in which case we require it to be independent of the shock sequence $ w_1, \\ldots, w_T $.  \n",
    "\n",
    "\n",
    "\n",
    "<a id='lq-cp'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4eda38db",
   "metadata": {},
   "source": [
    "### Information\n",
    "\n",
    "There’s one constraint we’ve neglected to mention so far, which is that the\n",
    "decision-maker who solves this LQ problem knows only the present and the past,\n",
    "not the future.\n",
    "\n",
    "To clarify this point, consider the sequence of controls $ \\{u_0, \\ldots, u_{T-1}\\} $.\n",
    "\n",
    "When choosing these controls, the decision-maker is permitted to take into account the effects of the shocks\n",
    "$ \\{w_1, \\ldots, w_T\\} $ on the system.\n",
    "\n",
    "However, it is typically assumed — and will be assumed here — that the\n",
    "time-$ t $ control $ u_t $ can  be made with knowledge of past and\n",
    "present shocks only.\n",
    "\n",
    "The fancy [measure-theoretic](https://en.wikipedia.org/wiki/Measure_%28mathematics%29) way of saying this is that $ u_t $ must be measurable with respect to the $ \\sigma $-algebra generated by $ x_0, w_1, w_2,\n",
    "\\ldots, w_t $.\n",
    "\n",
    "This is in fact equivalent to stating that $ u_t $ can be written in the form $ u_t = g_t(x_0, w_1, w_2, \\ldots, w_t) $ for some Borel measurable function $ g_t $.\n",
    "\n",
    "(Just about every function that’s useful for applications is Borel measurable,\n",
    "so, for the purposes of intuition, you can read that last phrase as “for some function $ g_t $”)\n",
    "\n",
    "Now note that $ x_t $ will ultimately depend on the realizations of $ x_0, w_1, w_2, \\ldots, w_t $.\n",
    "\n",
    "In fact, it turns out that $ x_t $ summarizes all the information about  these historical  shocks that the decision-maker needs to set controls optimally.\n",
    "\n",
    "More precisely, it can be shown that any optimal control $ u_t $ can always be written as a function of the current state alone.\n",
    "\n",
    "Hence in what follows we restrict attention to control policies (i.e., functions) of the form $ u_t = g_t(x_t) $.\n",
    "\n",
    "Actually, the preceding discussion applies to all standard dynamic programming problems.\n",
    "\n",
    "What’s special about the LQ case is that – as we shall soon see —  the optimal $ u_t $ turns out to be a linear function of $ x_t $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43e57659",
   "metadata": {},
   "source": [
    "### Solution\n",
    "\n",
    "To solve the finite horizon LQ problem we can use a dynamic programming\n",
    "strategy based on backward induction that is conceptually similar to the approach adopted in [this lecture](https://intro.quantecon.org/short_path.html).\n",
    "\n",
    "For reasons that will soon become clear, we first introduce the notation $ J_T(x) = x' R_f x $.\n",
    "\n",
    "Now consider the problem of the decision-maker in the second to last period.\n",
    "\n",
    "In particular, let the time be $ T-1 $, and suppose that the\n",
    "state is $ x_{T-1} $.\n",
    "\n",
    "The decision-maker must trade-off current and (discounted) final losses, and hence\n",
    "solves\n",
    "\n",
    "$$\n",
    "\\min_u \\{\n",
    "x_{T-1}' R x_{T-1} + u' Q u + \\beta \\,\n",
    "\\mathbb E J_T(A x_{T-1} + B u + C w_T)\n",
    "\\}\n",
    "$$\n",
    "\n",
    "At this stage, it is convenient to define the function\n",
    "\n",
    "\n",
    "<a id='equation-lq-lsm'></a>\n",
    "$$\n",
    "J_{T-1} (x) =\n",
    "\\min_u \\{\n",
    "x' R x + u' Q u + \\beta \\,\n",
    "\\mathbb E J_T(A x + B u + C w_T)\n",
    "\\} \\tag{56.7}\n",
    "$$\n",
    "\n",
    "The function $ J_{T-1} $ will be called the $ T-1 $ value function, and $ J_{T-1}(x) $ can be thought of as representing total “loss-to-go” from state $ x $ at time $ T-1 $ when the decision-maker behaves optimally.\n",
    "\n",
    "Now let’s step back to $ T-2 $.\n",
    "\n",
    "For a decision-maker at $ T-2 $, the value $ J_{T-1}(x) $ plays a role analogous to that played by the terminal loss $ J_T(x) = x' R_f x $ for the decision-maker at $ T-1 $.\n",
    "\n",
    "That is, $ J_{T-1}(x) $ summarizes the future loss associated with moving to state $ x $.\n",
    "\n",
    "The decision-maker chooses her control $ u $ to trade off current loss against future loss, where\n",
    "\n",
    "- the next period state is $ x_{T-1} = Ax_{T-2} + B u + C w_{T-1} $, and hence depends on the choice of current control.  \n",
    "- the “cost” of landing in state $ x_{T-1} $ is $ J_{T-1}(x_{T-1}) $.  \n",
    "\n",
    "\n",
    "Her problem is therefore\n",
    "\n",
    "$$\n",
    "\\min_u\n",
    "\\{\n",
    "x_{T-2}' R x_{T-2} + u' Q u + \\beta \\,\n",
    "\\mathbb E J_{T-1}(Ax_{T-2} + B u + C w_{T-1})\n",
    "\\}\n",
    "$$\n",
    "\n",
    "Letting\n",
    "\n",
    "$$\n",
    "J_{T-2} (x)\n",
    "= \\min_u\n",
    "\\{\n",
    "x' R x + u' Q u + \\beta \\,\n",
    "\\mathbb E J_{T-1}(Ax + B u + C w_{T-1})\n",
    "\\}\n",
    "$$\n",
    "\n",
    "the pattern for backward induction is now clear.\n",
    "\n",
    "In particular, we define a sequence of value functions $ \\{J_0, \\ldots, J_T\\} $ via\n",
    "\n",
    "$$\n",
    "J_{t-1} (x)\n",
    "= \\min_u\n",
    "\\{\n",
    "x' R x + u' Q u + \\beta \\,\n",
    "\\mathbb E J_{t}(Ax + B u + C w_t)\n",
    "\\}\n",
    "\\quad \\text{and} \\quad\n",
    "J_T(x) = x' R_f x\n",
    "$$\n",
    "\n",
    "The first equality is the Bellman equation from dynamic programming theory specialized to the finite horizon LQ problem.\n",
    "\n",
    "Now that we have $ \\{J_0, \\ldots, J_T\\} $, we can obtain the optimal controls.\n",
    "\n",
    "As a first step, let’s find out what the value functions look like.\n",
    "\n",
    "It turns out that every $ J_t $ has the form $ J_t(x) = x' P_t x + d_t $ where $ P_t $ is a $ n \\times n $ matrix and $ d_t $ is a constant.\n",
    "\n",
    "We can show this by induction, starting from $ P_T := R_f $ and $ d_T = 0 $.\n",
    "\n",
    "Using this notation, [(56.7)](#equation-lq-lsm) becomes\n",
    "\n",
    "\n",
    "<a id='equation-lq-fswb'></a>\n",
    "$$\n",
    "J_{T-1} (x) =\n",
    "\\min_u \\{\n",
    "x' R x + u' Q u + \\beta \\,\n",
    "\\mathbb E (A x + B u + C w_T)' P_T (A x + B u + C w_T)\n",
    "\\} \\tag{56.8}\n",
    "$$\n",
    "\n",
    "To obtain the minimizer, we can take the derivative of the r.h.s. with respect to $ u $ and set it equal to zero.\n",
    "\n",
    "Applying the relevant rules of [matrix calculus](https://python.quantecon.org/linear_algebra.html#la-mcalc), this gives\n",
    "\n",
    "\n",
    "<a id='equation-lq-oc0'></a>\n",
    "$$\n",
    "u  = - (Q + \\beta B' P_T B)^{-1} \\beta B' P_T A x \\tag{56.9}\n",
    "$$\n",
    "\n",
    "Plugging this back into [(56.8)](#equation-lq-fswb) and rearranging yields\n",
    "\n",
    "$$\n",
    "J_{T-1} (x) = x' P_{T-1} x + d_{T-1}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "\n",
    "<a id='equation-lq-finr'></a>\n",
    "$$\n",
    "P_{T-1} = R - \\beta^2 A' P_T B (Q + \\beta B' P_T B)^{-1} B' P_T A +\n",
    "\\beta A' P_T A \\tag{56.10}\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "\n",
    "<a id='equation-lq-finrd'></a>\n",
    "$$\n",
    "d_{T-1} := \\beta \\mathop{\\mathrm{trace}}(C' P_T C) \\tag{56.11}\n",
    "$$\n",
    "\n",
    "(The algebra is a good exercise — we’ll leave it up to you.)\n",
    "\n",
    "If we continue working backwards in this manner, it soon becomes clear that $ J_t (x) = x' P_t x + d_t $ as claimed, where $ \\{P_t\\} $ and $ \\{d_t\\} $ satisfy the recursions\n",
    "\n",
    "\n",
    "<a id='equation-lq-pr'></a>\n",
    "$$\n",
    "P_{t-1} = R - \\beta^2 A' P_t B (Q + \\beta B' P_t B)^{-1} B' P_t A +\n",
    "\\beta A' P_t A\n",
    "\\quad \\text{with } \\quad\n",
    "P_T = R_f \\tag{56.12}\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "\n",
    "<a id='equation-lq-dd'></a>\n",
    "$$\n",
    "d_{t-1} = \\beta (d_t + \\mathop{\\mathrm{trace}}(C' P_t C))\n",
    "\\quad \\text{with } \\quad\n",
    "d_T = 0 \\tag{56.13}\n",
    "$$\n",
    "\n",
    "Recalling [(56.9)](#equation-lq-oc0), the minimizers from these backward steps are\n",
    "\n",
    "\n",
    "<a id='equation-lq-oc'></a>\n",
    "$$\n",
    "u_t  = - F_t x_t\n",
    "\\quad \\text{where} \\quad\n",
    "F_t := (Q + \\beta B' P_{t+1} B)^{-1} \\beta B' P_{t+1} A \\tag{56.14}\n",
    "$$\n",
    "\n",
    "These are the linear optimal control policies we [discussed above](#lq-cp).\n",
    "\n",
    "In particular,  the sequence of controls given by [(56.14)](#equation-lq-oc) and [(56.1)](#equation-lq-lom) solves our finite horizon LQ problem.\n",
    "\n",
    "Rephrasing this more precisely, the sequence $ u_0, \\ldots, u_{T-1} $ given by\n",
    "\n",
    "\n",
    "<a id='equation-lq-xud'></a>\n",
    "$$\n",
    "u_t = - F_t x_t\n",
    "\\quad \\text{with} \\quad\n",
    "x_{t+1} = (A - BF_t) x_t + C w_{t+1} \\tag{56.15}\n",
    "$$\n",
    "\n",
    "for $ t = 0, \\ldots, T-1 $ attains the minimum of [(56.6)](#equation-lq-object) subject to our constraints."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b079a155",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "We will use code from [lqcontrol.py](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lqcontrol.py)\n",
    "in [QuantEcon.py](http://quantecon.org/quantecon-py)\n",
    "to solve finite and infinite horizon linear quadratic control problems.\n",
    "\n",
    "In the module, the various updating, simulation and fixed point methods\n",
    "are wrapped in a class  called `LQ`, which includes\n",
    "\n",
    "- Instance data:  \n",
    "  - The required parameters $ Q, R, A, B $ and optional parameters $ C, \\beta, T, R_f, N $ specifying a given LQ model  \n",
    "    - set $ T $ and $ R_f $ to `None` in the infinite horizon case  \n",
    "    - set `C = None` (or zero) in the deterministic case  \n",
    "  - the value function and policy data  \n",
    "    - $ d_t, P_t, F_t $ in the finite horizon case  \n",
    "    - $ d, P, F $ in the infinite horizon case  \n",
    "- Methods:  \n",
    "  - `update_values` — shifts $ d_t, P_t, F_t $ to their $ t-1 $ values via [(56.12)](#equation-lq-pr), [(56.13)](#equation-lq-dd) and [(56.14)](#equation-lq-oc)  \n",
    "  - `stationary_values` — computes $ P, d, F $ in the infinite horizon case  \n",
    "  - `compute_sequence` —- simulates the dynamics of $ x_t, u_t, w_t $ given $ x_0 $ and assuming standard normal shocks  \n",
    "\n",
    "\n",
    "\n",
    "<a id='lq-mfpa'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d63197a6",
   "metadata": {},
   "source": [
    "### An Application\n",
    "\n",
    "Early Keynesian models assumed that households have a constant marginal\n",
    "propensity to consume from current income.\n",
    "\n",
    "Data contradicted the constancy of the marginal propensity to consume.\n",
    "\n",
    "In response, Milton Friedman, Franco Modigliani and others built models\n",
    "based on a consumer’s preference for an intertemporally smooth consumption stream.\n",
    "\n",
    "(See, for example, [[Friedman, 1956](https://python.quantecon.org/zreferences.html#id164)] or [[Modigliani and Brumberg, 1954](https://python.quantecon.org/zreferences.html#id199)].)\n",
    "\n",
    "One property of those models is that households purchase and sell financial assets to make consumption streams smoother than income streams.\n",
    "\n",
    "The household savings problem [outlined above](#lq-hhp) captures these ideas.\n",
    "\n",
    "The optimization problem for the household is to choose a consumption sequence in order to minimize\n",
    "\n",
    "\n",
    "<a id='equation-lq-pio'></a>\n",
    "$$\n",
    "\\mathbb E \\,\n",
    "\\left\\{\n",
    "    \\sum_{t=0}^{T-1} \\beta^t (c_t - \\bar c)^2 + \\beta^T q a_T^2\n",
    "\\right\\} \\tag{56.16}\n",
    "$$\n",
    "\n",
    "subject to the sequence of budget constraints $ a_{t+1} = (1 + r) a_t - c_t + y_t, \\ t \\geq 0 $.\n",
    "\n",
    "Here $ q $ is a large positive constant, the role of which is to induce the consumer to target zero debt at the end of her life.\n",
    "\n",
    "(Without such a constraint, the optimal choice is to choose $ c_t = \\bar c $ in each period, letting assets adjust accordingly.)\n",
    "\n",
    "As before we set $ y_t = \\sigma w_{t+1} + \\mu $ and $ u_t := c_t - \\bar c $, after which the constraint can be written as in [(56.2)](#equation-lq-lomwc).\n",
    "\n",
    "We saw how this constraint could be manipulated into the LQ formulation $ x_{t+1} =\n",
    "Ax_t + Bu_t + Cw_{t+1} $ by setting $ x_t = (a_t \\; 1)' $ and using the definitions in [(56.4)](#equation-lq-lowmc2).\n",
    "\n",
    "To match with this state and control, the objective function [(56.16)](#equation-lq-pio) can\n",
    "be written in the form of [(56.6)](#equation-lq-object) by choosing\n",
    "\n",
    "$$\n",
    "Q := 1,\n",
    "\\quad\n",
    "R :=\n",
    "\\left(\n",
    "\\begin{array}{cc}\n",
    "0 & 0 \\\\\n",
    "0 & 0\n",
    "\\end{array}\n",
    "\\right),\n",
    "\\quad \\text{and} \\quad\n",
    "R_f :=\n",
    "\\left(\n",
    "\\begin{array}{cc}\n",
    "q & 0 \\\\\n",
    "0 & 0\n",
    "\\end{array}\n",
    "\\right)\n",
    "$$\n",
    "\n",
    "Now that the problem is expressed in LQ form, we can proceed to the solution\n",
    "by applying [(56.12)](#equation-lq-pr) and [(56.14)](#equation-lq-oc).\n",
    "\n",
    "After generating shocks $ w_1, \\ldots, w_T $, the dynamics for assets and\n",
    "consumption can be simulated via [(56.15)](#equation-lq-xud).\n",
    "\n",
    "The following figure was computed using $ r = 0.05, \\beta = 1 / (1+ r),\n",
    "\\bar c = 2,  \\mu = 1, \\sigma = 0.25, T = 45 $ and $ q = 10^6 $.\n",
    "\n",
    "The shocks $ \\{w_t\\} $ were taken to be IID and standard normal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4474126",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Model parameters\n",
    "r = 0.05\n",
    "β = 1/(1 + r)\n",
    "T = 45\n",
    "c_bar = 2\n",
    "σ = 0.25\n",
    "μ = 1\n",
    "q = 1e6\n",
    "\n",
    "# Formulate as an LQ problem\n",
    "Q = 1\n",
    "R = np.zeros((2, 2))\n",
    "Rf = np.zeros((2, 2))\n",
    "Rf[0, 0] = q\n",
    "A = [[1 + r, -c_bar + μ],\n",
    "    [0,              1]]\n",
    "B = [[-1],\n",
    "    [ 0]]\n",
    "C = [[σ],\n",
    "    [0]]\n",
    "\n",
    "# Compute solutions and simulate\n",
    "lq = LQ(Q, R, A, B, C, beta=β, T=T, Rf=Rf)\n",
    "x0 = (0, 1)\n",
    "xp, up, wp = lq.compute_sequence(x0)\n",
    "\n",
    "# Convert back to assets, consumption and income\n",
    "assets = xp[0, :]           # a_t\n",
    "c = up.flatten() + c_bar    # c_t\n",
    "income = σ * wp[0, 1:] + μ  # y_t\n",
    "\n",
    "# Plot results\n",
    "n_rows = 2\n",
    "fig, axes = plt.subplots(n_rows, 1, figsize=(12, 10))\n",
    "\n",
    "plt.subplots_adjust(hspace=0.5)\n",
    "\n",
    "bbox = (0., 1.02, 1., .102)\n",
    "legend_args = {'bbox_to_anchor': bbox, 'loc': 3, 'mode': 'expand'}\n",
    "p_args = {'lw': 2, 'alpha': 0.7}\n",
    "\n",
    "axes[0].plot(list(range(1, T+1)), income, 'g-', label=\"non-financial income\",\n",
    "            **p_args)\n",
    "axes[0].plot(list(range(T)), c, 'k-', label=\"consumption\", **p_args)\n",
    "\n",
    "axes[1].plot(list(range(1, T+1)), np.cumsum(income - μ), 'r-',\n",
    "            label=\"cumulative unanticipated income\", **p_args)\n",
    "axes[1].plot(list(range(T+1)), assets, 'b-', label=\"assets\", **p_args)\n",
    "axes[1].plot(list(range(T)), np.zeros(T), 'k-')\n",
    "\n",
    "for ax in axes:\n",
    "    ax.grid()\n",
    "    ax.set_xlabel('Time')\n",
    "    ax.legend(ncol=2, **legend_args)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6297f0f",
   "metadata": {},
   "source": [
    "The top panel shows the time path of consumption $ c_t $ and income $ y_t $ in the simulation.\n",
    "\n",
    "As anticipated by the discussion on consumption smoothing, the time path of\n",
    "consumption is much smoother than that for income.\n",
    "\n",
    "(But note that  consumption becomes more irregular towards the end of life,\n",
    "when the zero final asset requirement impinges more on consumption choices.)\n",
    "\n",
    "The second panel in the figure shows that the time path of assets $ a_t $ is\n",
    "closely correlated with cumulative unanticipated income, where the latter is defined as\n",
    "\n",
    "$$\n",
    "z_t := \\sum_{j=0}^t \\sigma w_t\n",
    "$$\n",
    "\n",
    "A key message is that unanticipated windfall gains are saved rather\n",
    "than consumed, while unanticipated negative shocks are met by reducing assets.\n",
    "\n",
    "(Again, this relationship breaks down towards the end of life due to the zero final asset requirement.)\n",
    "\n",
    "These results are relatively robust to changes in parameters.\n",
    "\n",
    "For example, let’s increase $ \\beta $ from $ 1 / (1 + r) \\approx 0.952 $ to $ 0.96 $ while keeping other parameters fixed.\n",
    "\n",
    "This consumer is slightly more patient than the last one, and hence puts\n",
    "relatively more weight on later consumption values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee80eabb",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Compute solutions and simulate\n",
    "lq = LQ(Q, R, A, B, C, beta=0.96, T=T, Rf=Rf)\n",
    "x0 = (0, 1)\n",
    "xp, up, wp = lq.compute_sequence(x0)\n",
    "\n",
    "# Convert back to assets, consumption and income\n",
    "assets = xp[0, :]           # a_t\n",
    "c = up.flatten() + c_bar    # c_t\n",
    "income = σ * wp[0, 1:] + μ  # y_t\n",
    "\n",
    "# Plot results\n",
    "n_rows = 2\n",
    "fig, axes = plt.subplots(n_rows, 1, figsize=(12, 10))\n",
    "\n",
    "plt.subplots_adjust(hspace=0.5)\n",
    "\n",
    "bbox = (0., 1.02, 1., .102)\n",
    "legend_args = {'bbox_to_anchor': bbox, 'loc': 3, 'mode': 'expand'}\n",
    "p_args = {'lw': 2, 'alpha': 0.7}\n",
    "\n",
    "axes[0].plot(list(range(1, T+1)), income, 'g-', label=\"non-financial income\",\n",
    "             **p_args)\n",
    "axes[0].plot(list(range(T)), c, 'k-', label=\"consumption\", **p_args)\n",
    "\n",
    "axes[1].plot(list(range(1, T+1)), np.cumsum(income - μ), 'r-',\n",
    "             label=\"cumulative unanticipated income\", **p_args)\n",
    "axes[1].plot(list(range(T+1)), assets, 'b-', label=\"assets\", **p_args)\n",
    "axes[1].plot(list(range(T)), np.zeros(T), 'k-')\n",
    "\n",
    "for ax in axes:\n",
    "    ax.grid()\n",
    "    ax.set_xlabel('Time')\n",
    "    ax.legend(ncol=2, **legend_args)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af49e13b",
   "metadata": {},
   "source": [
    "We now have a slowly rising consumption stream and a hump-shaped build-up\n",
    "of assets in the middle periods to fund rising consumption.\n",
    "\n",
    "However, the essential features are the same: consumption is smooth relative to income, and assets are strongly positively correlated with cumulative unanticipated income."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6360bc05",
   "metadata": {},
   "source": [
    "## Extensions and Comments\n",
    "\n",
    "Let’s now consider a number of standard extensions to the LQ problem treated above."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca633dbc",
   "metadata": {},
   "source": [
    "### Time-Varying Parameters\n",
    "\n",
    "In some settings, it can be desirable to allow $ A, B, C, R $ and $ Q $ to depend on $ t $.\n",
    "\n",
    "For the sake of simplicity, we’ve chosen not to treat this extension in our implementation given below.\n",
    "\n",
    "However, the loss of generality is not as large as you might first imagine.\n",
    "\n",
    "In fact, we can tackle many models with time-varying parameters by suitable choice of state variables.\n",
    "\n",
    "One illustration is given [below](#lq-nsi).\n",
    "\n",
    "For further examples and a more systematic treatment, see [[Hansen and Sargent, 2013](https://python.quantecon.org/zreferences.html#id169)], section 2.4.\n",
    "\n",
    "\n",
    "<a id='lq-cpt'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "561baeb9",
   "metadata": {},
   "source": [
    "### Adding a Cross-Product Term\n",
    "\n",
    "In some LQ problems, preferences include a cross-product term $ u_t' N x_t $, so that the objective function becomes\n",
    "\n",
    "\n",
    "<a id='equation-lq-object-cp'></a>\n",
    "$$\n",
    "\\mathbb E \\,\n",
    "\\left\\{\n",
    "    \\sum_{t=0}^{T-1} \\beta^t (x_t' R x_t + u_t' Q u_t + 2 u_t' N x_t) + \\beta^T x_T' R_f x_T\n",
    "\\right\\} \\tag{56.17}\n",
    "$$\n",
    "\n",
    "Our results extend to this case in a straightforward way.\n",
    "\n",
    "The sequence $ \\{P_t\\} $ from [(56.12)](#equation-lq-pr) becomes\n",
    "\n",
    "\n",
    "<a id='equation-lq-pr-cp'></a>\n",
    "$$\n",
    "P_{t-1} = R - (\\beta B' P_t A + N)'\n",
    "(Q + \\beta B' P_t B)^{-1} (\\beta B' P_t A + N) +\n",
    "\\beta A' P_t A\n",
    "\\quad \\text{with } \\quad\n",
    "P_T = R_f \\tag{56.18}\n",
    "$$\n",
    "\n",
    "The policies in [(56.14)](#equation-lq-oc) are modified to\n",
    "\n",
    "\n",
    "<a id='equation-lq-oc-cp'></a>\n",
    "$$\n",
    "u_t  = - F_t x_t\n",
    "\\quad \\text{where} \\quad\n",
    "F_t := (Q + \\beta B' P_{t+1} B)^{-1} (\\beta B' P_{t+1} A + N) \\tag{56.19}\n",
    "$$\n",
    "\n",
    "The sequence $ \\{d_t\\} $ is unchanged from [(56.13)](#equation-lq-dd).\n",
    "\n",
    "We leave interested readers to confirm these results (the calculations are long but not overly difficult).\n",
    "\n",
    "\n",
    "<a id='lq-ih'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "036a5efc",
   "metadata": {},
   "source": [
    "### Infinite Horizon\n",
    "\n",
    "\n",
    "<a id='index-2'></a>\n",
    "Finally, we consider the infinite horizon case, with [cross-product term](#lq-cpt), unchanged dynamics and\n",
    "objective function given by\n",
    "\n",
    "\n",
    "<a id='equation-lq-object-ih'></a>\n",
    "$$\n",
    "\\mathbb E \\,\n",
    "\\left\\{\n",
    "    \\sum_{t=0}^{\\infty} \\beta^t (x_t' R x_t + u_t' Q u_t + 2 u_t' N x_t)\n",
    "\\right\\} \\tag{56.20}\n",
    "$$\n",
    "\n",
    "In the infinite horizon case, optimal policies can depend on time\n",
    "only if time itself is a component of the  state vector $ x_t $.\n",
    "\n",
    "In other words, there exists a fixed matrix $ F $ such that $ u_t = -\n",
    "F x_t $ for all $ t $.\n",
    "\n",
    "That decision rules are constant over time is intuitive — after all, the decision-maker faces the\n",
    "same infinite horizon at every stage, with only the current state changing.\n",
    "\n",
    "Not surprisingly, $ P $ and $ d $ are also constant.\n",
    "\n",
    "The stationary matrix $ P $ is the solution to the\n",
    "[discrete-time algebraic Riccati equation](https://en.wikipedia.org/wiki/Algebraic_Riccati_equation).\n",
    "\n",
    "\n",
    "<a id='riccati-equation'></a>\n",
    "\n",
    "<a id='equation-lq-pr-ih'></a>\n",
    "$$\n",
    "P = R - (\\beta B' P A + N)'\n",
    "(Q + \\beta B' P B)^{-1} (\\beta B' P A + N) +\n",
    "\\beta A' P A \\tag{56.21}\n",
    "$$\n",
    "\n",
    "Equation [(56.21)](#equation-lq-pr-ih) is also called the *LQ Bellman equation*, and the map\n",
    "that sends a given $ P $ into the right-hand side of [(56.21)](#equation-lq-pr-ih) is\n",
    "called the *LQ Bellman operator*.\n",
    "\n",
    "The stationary optimal policy for this model is\n",
    "\n",
    "\n",
    "<a id='equation-lq-oc-ih'></a>\n",
    "$$\n",
    "u  = - F x\n",
    "\\quad \\text{where} \\quad\n",
    "F = (Q + \\beta B' P B)^{-1} (\\beta B' P A + N) \\tag{56.22}\n",
    "$$\n",
    "\n",
    "The sequence $ \\{d_t\\} $ from [(56.13)](#equation-lq-dd) is replaced by the constant value\n",
    "\n",
    "\n",
    "<a id='equation-lq-dd-ih'></a>\n",
    "$$\n",
    "d\n",
    ":= \\mathop{\\mathrm{trace}}(C' P C) \\frac{\\beta}{1 - \\beta} \\tag{56.23}\n",
    "$$\n",
    "\n",
    "The state evolves according to the time-homogeneous process $ x_{t+1} = (A - BF) x_t + C w_{t+1} $.\n",
    "\n",
    "An example infinite horizon problem is treated [below](#lqc-mwac).\n",
    "\n",
    "\n",
    "<a id='lq-cert-eq'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56c074e5",
   "metadata": {},
   "source": [
    "### Certainty Equivalence\n",
    "\n",
    "Linear quadratic control problems of the class discussed above have the property of *certainty equivalence*.\n",
    "\n",
    "By this, we mean that the optimal policy $ F $ is not affected by the parameters in $ C $, which specify the shock process.\n",
    "\n",
    "This can be confirmed by inspecting [(56.22)](#equation-lq-oc-ih) or [(56.19)](#equation-lq-oc-cp).\n",
    "\n",
    "It follows that we can ignore uncertainty when solving for optimal behavior, and plug it back in when examining optimal state dynamics."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dec0b633",
   "metadata": {},
   "source": [
    "## Further Applications\n",
    "\n",
    "\n",
    "<a id='lq-nsi'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb695f64",
   "metadata": {},
   "source": [
    "### Application 1: Age-Dependent Income Process\n",
    "\n",
    "[Previously](#lq-mfpa) we studied a permanent income model that generated consumption smoothing.\n",
    "\n",
    "One unrealistic feature of that model is the assumption that the mean of the random income process does not depend on the consumer’s age.\n",
    "\n",
    "A more realistic income profile is one that rises in early working life, peaks towards the middle and maybe declines toward the end of working life and falls more during retirement.\n",
    "\n",
    "In this section, we will model this rise and fall as a symmetric inverted “U” using a polynomial in age.\n",
    "\n",
    "As before, the consumer seeks to minimize\n",
    "\n",
    "\n",
    "<a id='equation-lq-pip'></a>\n",
    "$$\n",
    "\\mathbb E \\,\n",
    "\\left\\{\n",
    "    \\sum_{t=0}^{T-1} \\beta^t (c_t - \\bar c)^2 + \\beta^T q a_T^2\n",
    "\\right\\} \\tag{56.24}\n",
    "$$\n",
    "\n",
    "subject to $ a_{t+1} = (1 + r) a_t - c_t + y_t, \\ t \\geq 0 $.\n",
    "\n",
    "For income we now take $ y_t = p(t) + \\sigma w_{t+1} $ where $ p(t) := m_0 + m_1 t + m_2 t^2 $.\n",
    "\n",
    "(In [the next section](#lq-nsi2) we employ some tricks to implement a more sophisticated model.)\n",
    "\n",
    "The coefficients $ m_0, m_1, m_2 $ are chosen such that $ p(0)=0, p(T/2) = \\mu, $ and $ p(T)=0 $.\n",
    "\n",
    "You can confirm that the specification $ m_0 = 0, m_1 = T \\mu / (T/2)^2, m_2 = - \\mu / (T/2)^2 $ satisfies these constraints.\n",
    "\n",
    "To put this into an LQ setting, consider the budget constraint, which becomes\n",
    "\n",
    "\n",
    "<a id='equation-lq-hib'></a>\n",
    "$$\n",
    "a_{t+1} = (1 + r) a_t - u_t - \\bar c + m_1 t + m_2 t^2 + \\sigma w_{t+1} \\tag{56.25}\n",
    "$$\n",
    "\n",
    "The fact that $ a_{t+1} $ is a linear function of\n",
    "$ (a_t, 1, t, t^2) $ suggests taking these four variables as the state\n",
    "vector $ x_t $.\n",
    "\n",
    "Once a good choice of state and control (recall $ u_t = c_t - \\bar c $)\n",
    "has been made, the remaining specifications fall into place relatively easily.\n",
    "\n",
    "Thus, for the dynamics we set\n",
    "\n",
    "\n",
    "<a id='equation-lq-lowmc3'></a>\n",
    "$$\n",
    "x_t :=\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "a_t \\\\\n",
    "1 \\\\\n",
    "t \\\\\n",
    "t^2\n",
    "\\end{array}\n",
    "\\right),\n",
    "\\quad\n",
    "A :=\n",
    "\\left(\n",
    "\\begin{array}{cccc}\n",
    "1 + r & -\\bar c & m_1 & m_2 \\\\\n",
    "0     & 1       & 0   & 0   \\\\\n",
    "0     & 1       & 1   & 0   \\\\\n",
    "0     & 1       & 2   & 1\n",
    "\\end{array}\n",
    "\\right),\n",
    "\\quad\n",
    "B :=\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "-1 \\\\\n",
    "0 \\\\\n",
    "0 \\\\\n",
    "0\n",
    "\\end{array}\n",
    "\\right),\n",
    "\\quad\n",
    "C :=\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "\\sigma \\\\\n",
    "0 \\\\\n",
    "0 \\\\\n",
    "0\n",
    "\\end{array}\n",
    "\\right) \\tag{56.26}\n",
    "$$\n",
    "\n",
    "If you expand the expression $ x_{t+1} = A x_t + B u_t + C w_{t+1} $ using\n",
    "this specification, you will find that assets follow [(56.25)](#equation-lq-hib) as desired\n",
    "and that the other state variables also update appropriately.\n",
    "\n",
    "To implement preference specification [(56.24)](#equation-lq-pip) we take\n",
    "\n",
    "\n",
    "<a id='equation-lq-4sp'></a>\n",
    "$$\n",
    "Q := 1,\n",
    "\\quad\n",
    "R :=\n",
    "\\left(\n",
    "\\begin{array}{cccc}\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\quad \\text{and} \\quad\n",
    "R_f :=\n",
    "\\left(\n",
    "\\begin{array}{cccc}\n",
    "q & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0\n",
    "\\end{array}\n",
    "\\right) \\tag{56.27}\n",
    "$$\n",
    "\n",
    "The next figure shows a simulation of consumption and assets computed using\n",
    "the `compute_sequence` method of `lqcontrol.py` with initial assets set to zero.\n",
    "\n",
    "\n",
    "<a id='solution-lqc-ex1-fig'></a>\n",
    "![https://python.quantecon.org/_static/lecture_specific/lqcontrol/solution_lqc_ex1.png](https://python.quantecon.org/_static/lecture_specific/lqcontrol/solution_lqc_ex1.png)\n",
    "\n",
    "  \n",
    "Once again, smooth consumption is a dominant feature of the sample  paths.\n",
    "\n",
    "The asset path exhibits dynamics consistent with standard life cycle theory.\n",
    "\n",
    "Exercise 56.1 gives the full set of parameters used here and asks you to replicate the figure.\n",
    "\n",
    "\n",
    "<a id='lq-nsi2'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36971a02",
   "metadata": {},
   "source": [
    "### Application 2: A Permanent Income Model with Retirement\n",
    "\n",
    "In the [previous application](#lq-nsi), we generated income dynamics with an inverted U shape using polynomials and placed them in an LQ framework.\n",
    "\n",
    "It is arguably the case that this income process still contains unrealistic features.\n",
    "\n",
    "A more common earning profile is where\n",
    "\n",
    "1. income grows over working life, fluctuating around an increasing trend, with growth flattening off in later years  \n",
    "1. retirement follows, with lower but relatively stable (non-financial) income  \n",
    "\n",
    "\n",
    "Letting $ K $ be the retirement date, we can express these income dynamics\n",
    "by\n",
    "\n",
    "\n",
    "<a id='equation-lq-cases'></a>\n",
    "$$\n",
    "y_t =\n",
    "\\begin{cases}\n",
    "p(t) + \\sigma w_{t+1} & \\quad \\text{if } t \\leq K  \\\\\n",
    "s                     & \\quad \\text{otherwise }\n",
    "\\end{cases} \\tag{56.28}\n",
    "$$\n",
    "\n",
    "Here\n",
    "\n",
    "- $ p(t) := m_1 t + m_2 t^2 $ with the coefficients $ m_1, m_2 $ chosen such that $ p(K) = \\mu $ and $ p(0) = p(2K)=0 $  \n",
    "- $ s $ is retirement income  \n",
    "\n",
    "\n",
    "We suppose that preferences are unchanged and given by [(56.16)](#equation-lq-pio).\n",
    "\n",
    "The budget constraint is also unchanged and given by $ a_{t+1} = (1 + r) a_t - c_t + y_t $.\n",
    "\n",
    "Our aim is to solve this problem and simulate paths using the LQ techniques described in this lecture.\n",
    "\n",
    "In fact, this is a nontrivial problem, as the kink in the dynamics [(56.28)](#equation-lq-cases) at $ K $ makes it very difficult to express the law of motion as a fixed-coefficient linear system.\n",
    "\n",
    "However, we can still use our LQ methods here by suitably linking two-component LQ problems.\n",
    "\n",
    "These two LQ problems describe the consumer’s behavior during her working life (`lq_working`) and retirement (`lq_retired`).\n",
    "\n",
    "(This is possible because, in the two separate periods of life, the respective income processes\n",
    "[polynomial trend and constant] each fit the LQ framework.)\n",
    "\n",
    "The basic idea is that although the whole problem is not a single time-invariant LQ problem, it is\n",
    "still a dynamic programming problem, and hence we can use appropriate Bellman equations at\n",
    "every stage.\n",
    "\n",
    "Based on this logic, we can\n",
    "\n",
    "1. solve `lq_retired` by the usual backward induction procedure, iterating back to the start of retirement.  \n",
    "1. take the start-of-retirement value function generated by this process, and use it as  the terminal condition $ R_f $ to feed into the `lq_working` specification.  \n",
    "1. solve `lq_working` by backward induction from this choice of $ R_f $, iterating back to the start of working life.  \n",
    "\n",
    "\n",
    "This process gives the entire life-time sequence of value functions and optimal policies.\n",
    "\n",
    "The next figure shows one simulation based on this procedure.\n",
    "\n",
    "\n",
    "<a id='solution-lqc-ex2-fig'></a>\n",
    "![https://python.quantecon.org/_static/lecture_specific/lqcontrol/solution_lqc_ex2.png](https://python.quantecon.org/_static/lecture_specific/lqcontrol/solution_lqc_ex2.png)\n",
    "\n",
    "  \n",
    "The full set of parameters used in the simulation is discussed in Exercise 56.2, where you are asked to replicate the figure.\n",
    "\n",
    "Once again, the dominant feature observable in the simulation is consumption\n",
    "smoothing.\n",
    "\n",
    "The asset path fits well with standard life cycle theory, with dissaving early\n",
    "in life followed by later saving.\n",
    "\n",
    "Assets peak at retirement and subsequently decline.\n",
    "\n",
    "\n",
    "<a id='lqc-mwac'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b28dce20",
   "metadata": {},
   "source": [
    "### Application 3: Monopoly with Adjustment Costs\n",
    "\n",
    "Consider a monopolist facing stochastic inverse demand function\n",
    "\n",
    "$$\n",
    "p_t = a_0 - a_1 q_t + d_t\n",
    "$$\n",
    "\n",
    "Here $ q_t $ is output, and the demand shock $ d_t $ follows\n",
    "\n",
    "$$\n",
    "d_{t+1} = \\rho d_t + \\sigma w_{t+1}\n",
    "$$\n",
    "\n",
    "where $ \\{w_t\\} $ is IID and standard normal.\n",
    "\n",
    "The monopolist maximizes the expected discounted sum of present and future profits\n",
    "\n",
    "\n",
    "<a id='equation-lq-object-mp'></a>\n",
    "$$\n",
    "\\mathbb E \\,\n",
    "\\left\\{\n",
    "    \\sum_{t=0}^{\\infty} \\beta^t\n",
    "    \\pi_t\n",
    "\\right\\}\n",
    "\\quad \\text{where} \\quad\n",
    "\\pi_t := p_t q_t - c q_t - \\gamma (q_{t+1} - q_t)^2 \\tag{56.29}\n",
    "$$\n",
    "\n",
    "Here\n",
    "\n",
    "- $ \\gamma (q_{t+1} - q_t)^2 $ represents adjustment costs  \n",
    "- $ c $ is average cost of production  \n",
    "\n",
    "\n",
    "This can be formulated as an LQ problem and then solved and simulated,\n",
    "but first let’s study the problem and try to get some intuition.\n",
    "\n",
    "One way to start thinking about the problem is to consider what would happen\n",
    "if $ \\gamma = 0 $.\n",
    "\n",
    "Without adjustment costs there is no intertemporal trade-off, so the\n",
    "monopolist will choose output to maximize current profit in each period.\n",
    "\n",
    "It’s not difficult to show that profit-maximizing output is\n",
    "\n",
    "$$\n",
    "\\bar q_t := \\frac{a_0 - c + d_t}{2 a_1}\n",
    "$$\n",
    "\n",
    "In light of this discussion, what we might expect for general $ \\gamma $ is that\n",
    "\n",
    "- if $ \\gamma $ is close to zero, then $ q_t $ will track the time path of $ \\bar q_t $ relatively closely.  \n",
    "- if $ \\gamma $ is larger, then $ q_t $ will be smoother than $ \\bar q_t $, as the monopolist seeks to avoid adjustment costs.  \n",
    "\n",
    "\n",
    "This intuition turns out to be correct.\n",
    "\n",
    "The following figures show simulations produced by solving the corresponding LQ problem.\n",
    "\n",
    "The only difference in parameters across the figures is the size of $ \\gamma $\n",
    "\n",
    "![https://python.quantecon.org/_static/lecture_specific/lqcontrol/solution_lqc_ex3_g1.png](https://python.quantecon.org/_static/lecture_specific/lqcontrol/solution_lqc_ex3_g1.png)\n",
    "\n",
    "  \n",
    "![https://python.quantecon.org/_static/lecture_specific/lqcontrol/solution_lqc_ex3_g10.png](https://python.quantecon.org/_static/lecture_specific/lqcontrol/solution_lqc_ex3_g10.png)\n",
    "\n",
    "  \n",
    "![https://python.quantecon.org/_static/lecture_specific/lqcontrol/solution_lqc_ex3_g50.png](https://python.quantecon.org/_static/lecture_specific/lqcontrol/solution_lqc_ex3_g50.png)\n",
    "\n",
    "  \n",
    "To produce these figures we converted the monopolist problem into an LQ problem.\n",
    "\n",
    "The key to this conversion is to choose the right state — which can be a bit of an art.\n",
    "\n",
    "Here we take $ x_t = (\\bar q_t \\;\\, q_t \\;\\, 1)' $, while the control is chosen as $ u_t = q_{t+1} - q_t $.\n",
    "\n",
    "We also manipulated the profit function slightly.\n",
    "\n",
    "In [(56.29)](#equation-lq-object-mp), current profits are $ \\pi_t := p_t q_t - c q_t - \\gamma (q_{t+1} - q_t)^2 $.\n",
    "\n",
    "Let’s now replace $ \\pi_t $ in [(56.29)](#equation-lq-object-mp) with $ \\hat \\pi_t := \\pi_t - a_1 \\bar q_t^2 $.\n",
    "\n",
    "This makes no difference to the solution, since $ a_1 \\bar q_t^2 $ does not depend on the controls.\n",
    "\n",
    "(In fact, we are just adding a constant term to [(56.29)](#equation-lq-object-mp), and optimizers are not affected by constant terms.)\n",
    "\n",
    "The reason for making this substitution is that, as you will be able to\n",
    "verify, $ \\hat \\pi_t $ reduces to the simple quadratic\n",
    "\n",
    "$$\n",
    "\\hat \\pi_t = -a_1 (q_t - \\bar q_t)^2 - \\gamma u_t^2\n",
    "$$\n",
    "\n",
    "After negation to convert to a minimization problem, the objective becomes\n",
    "\n",
    "\n",
    "<a id='equation-lq-object-mp2'></a>\n",
    "$$\n",
    "\\min\n",
    "\\mathbb E \\,\n",
    "    \\sum_{t=0}^{\\infty} \\beta^t\n",
    "\\left\\{\n",
    "    a_1 ( q_t - \\bar q_t)^2 + \\gamma u_t^2\n",
    "\\right\\} \\tag{56.30}\n",
    "$$\n",
    "\n",
    "It’s now relatively straightforward to find $ R $ and $ Q $ such that\n",
    "[(56.30)](#equation-lq-object-mp2) can be written as [(56.20)](#equation-lq-object-ih).\n",
    "\n",
    "Furthermore, the matrices $ A, B $ and $ C $ from [(56.1)](#equation-lq-lom)\n",
    "can be found by writing down the dynamics of each element of the state.\n",
    "\n",
    "Exercise 56.3 asks you to complete this process, and reproduce the preceding figures."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "614b3db7",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "413a8bd7",
   "metadata": {},
   "source": [
    "## Exercise 56.1\n",
    "\n",
    "Replicate the figure with polynomial income [shown above](#solution-lqc-ex1-fig).\n",
    "\n",
    "The parameters are $ r = 0.05, \\beta = 1 / (1 + r), \\bar c = 1.5,  \\mu = 2, \\sigma = 0.15, T = 50 $ and $ q = 10^4 $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6abe0f67",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 56.1](https://python.quantecon.org/#lqc_ex1)\n",
    "\n",
    "Here’s one solution.\n",
    "\n",
    "We use some fancy plot commands to get a certain style — feel free to\n",
    "use simpler ones.\n",
    "\n",
    "The model is an LQ permanent income / life-cycle model with hump-shaped\n",
    "income\n",
    "\n",
    "$$\n",
    "y_t = m_1 t + m_2 t^2 + \\sigma w_{t+1}\n",
    "$$\n",
    "\n",
    "where $ \\{w_t\\} $ is IID $ N(0, 1) $ and the coefficients\n",
    "$ m_1 $ and $ m_2 $ are chosen so that\n",
    "$ p(t) = m_1 t + m_2 t^2 $ has an inverted U shape with\n",
    "\n",
    "- $ p(0) = 0, p(T/2) = \\mu $, and  \n",
    "- $ p(T) = 0 $  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eebc7857",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Model parameters\n",
    "r = 0.05\n",
    "β = 1/(1 + r)\n",
    "T = 50\n",
    "c_bar = 1.5\n",
    "σ = 0.15\n",
    "μ = 2\n",
    "q = 1e4\n",
    "m1 = T * (μ/(T/2)**2)\n",
    "m2 = -(μ/(T/2)**2)\n",
    "\n",
    "# Formulate as an LQ problem\n",
    "Q = 1\n",
    "R = np.zeros((4, 4))\n",
    "Rf = np.zeros((4, 4))\n",
    "Rf[0, 0] = q\n",
    "A = [[1 + r, -c_bar, m1, m2],\n",
    "     [0,          1,  0,  0],\n",
    "     [0,          1,  1,  0],\n",
    "     [0,          1,  2,  1]]\n",
    "B = [[-1],\n",
    "     [ 0],\n",
    "     [ 0],\n",
    "     [ 0]]\n",
    "C = [[σ],\n",
    "     [0],\n",
    "     [0],\n",
    "     [0]]\n",
    "\n",
    "# Compute solutions and simulate\n",
    "lq = LQ(Q, R, A, B, C, beta=β, T=T, Rf=Rf)\n",
    "x0 = (0, 1, 0, 0)\n",
    "xp, up, wp = lq.compute_sequence(x0)\n",
    "\n",
    "# Convert results back to assets, consumption and income\n",
    "ap = xp[0, :]               # Assets\n",
    "c = up.flatten() + c_bar    # Consumption\n",
    "time = np.arange(1, T+1)\n",
    "income = σ * wp[0, 1:] + m1 * time + m2 * time**2  # Income\n",
    "\n",
    "\n",
    "# Plot results\n",
    "n_rows = 2\n",
    "fig, axes = plt.subplots(n_rows, 1, figsize=(12, 10))\n",
    "\n",
    "plt.subplots_adjust(hspace=0.5)\n",
    "\n",
    "bbox = (0., 1.02, 1., .102)\n",
    "legend_args = {'bbox_to_anchor': bbox, 'loc': 3, 'mode': 'expand'}\n",
    "p_args = {'lw': 2, 'alpha': 0.7}\n",
    "\n",
    "axes[0].plot(range(1, T+1), income, 'g-', label=\"non-financial income\",\n",
    "            **p_args)\n",
    "axes[0].plot(range(T), c, 'k-', label=\"consumption\", **p_args)\n",
    "\n",
    "axes[1].plot(range(T+1), ap.flatten(), 'b-', label=\"assets\", **p_args)\n",
    "axes[1].plot(range(T+1), np.zeros(T+1), 'k-')\n",
    "\n",
    "for ax in axes:\n",
    "    ax.grid()\n",
    "    ax.set_xlabel('Time')\n",
    "    ax.legend(ncol=2, **legend_args)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f7b7041",
   "metadata": {},
   "source": [
    "## Exercise 56.2\n",
    "\n",
    "Replicate the figure on work and retirement [shown above](#solution-lqc-ex2-fig).\n",
    "\n",
    "The parameters are $ r = 0.05, \\beta = 1 / (1 + r), \\bar c = 4,  \\mu = 4, \\sigma = 0.35, K = 40, T = 60, s = 1 $ and $ q = 10^4 $.\n",
    "\n",
    "To understand the overall procedure, carefully read the section containing that figure.\n",
    "\n",
    "First, in order to make our approach work, we must ensure that both LQ problems have the same state variables and control.\n",
    "\n",
    "As with previous applications, the control can be set to $ u_t = c_t - \\bar c $.\n",
    "\n",
    "For `lq_working`, $ x_t, A, B, C $ can be chosen as in [(56.26)](#equation-lq-lowmc3).\n",
    "\n",
    "- Recall that $ m_1, m_2 $ are chosen so that $ p(K) = \\mu $ and $ p(2K)=0 $.  \n",
    "\n",
    "\n",
    "For `lq_retired`, use the same definition of $ x_t $ and $ u_t $, but modify $ A, B, C $ to correspond to constant income $ y_t = s $.\n",
    "\n",
    "For `lq_retired`, set preferences as in [(56.27)](#equation-lq-4sp).\n",
    "\n",
    "For `lq_working`, preferences are the same, except that $ R_f $ should\n",
    "be replaced by the final value function that emerges from iterating `lq_retired`\n",
    "back to the start of retirement.\n",
    "\n",
    "With some careful footwork, the simulation can be generated by patching\n",
    "together the simulations from these two separate models."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eec712fd",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 56.2](https://python.quantecon.org/#lqc_ex2)\n",
    "\n",
    "This is a permanent income / life-cycle model with polynomial growth in\n",
    "income over working life followed by a fixed retirement income.\n",
    "\n",
    "The model is solved by combining two LQ programming problems as described in\n",
    "the lecture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fc64c29",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Model parameters\n",
    "r = 0.05\n",
    "β = 1/(1 + r)\n",
    "T = 60\n",
    "K = 40\n",
    "c_bar = 4\n",
    "σ = 0.35\n",
    "μ = 4\n",
    "q = 1e4\n",
    "s = 1\n",
    "m1 = 2 * μ/K\n",
    "m2 = -μ/K**2\n",
    "\n",
    "# Formulate LQ problem 1 (retirement)\n",
    "Q = 1\n",
    "R = np.zeros((4, 4))\n",
    "Rf = np.zeros((4, 4))\n",
    "Rf[0, 0] = q\n",
    "A = [[1 + r, s - c_bar, 0, 0],\n",
    "     [0,             1, 0, 0],\n",
    "     [0,             1, 1, 0],\n",
    "     [0,             1, 2, 1]]\n",
    "B = [[-1],\n",
    "     [ 0],\n",
    "     [ 0],\n",
    "     [ 0]]\n",
    "C = [[0],\n",
    "     [0],\n",
    "     [0],\n",
    "     [0]]\n",
    "\n",
    "# Initialize LQ instance for retired agent\n",
    "lq_retired = LQ(Q, R, A, B, C, beta=β, T=T-K, Rf=Rf)\n",
    "# Iterate back to start of retirement, record final value function\n",
    "for i in range(T-K):\n",
    "    lq_retired.update_values()\n",
    "Rf2 = lq_retired.P\n",
    "\n",
    "# Formulate LQ problem 2 (working life)\n",
    "R = np.zeros((4, 4))\n",
    "A = [[1 + r, -c_bar, m1, m2],\n",
    "     [0,          1,  0,  0],\n",
    "     [0,          1,  1,  0],\n",
    "     [0,          1,  2,  1]]\n",
    "B = [[-1],\n",
    "     [ 0],\n",
    "     [ 0],\n",
    "     [ 0]]\n",
    "C = [[σ],\n",
    "     [0],\n",
    "     [0],\n",
    "     [0]]\n",
    "\n",
    "# Set up working life LQ instance with terminal Rf from lq_retired\n",
    "lq_working = LQ(Q, R, A, B, C, beta=β, T=K, Rf=Rf2)\n",
    "\n",
    "# Simulate working state / control paths\n",
    "x0 = (0, 1, 0, 0)\n",
    "xp_w, up_w, wp_w = lq_working.compute_sequence(x0)\n",
    "# Simulate retirement paths (note the initial condition)\n",
    "xp_r, up_r, wp_r = lq_retired.compute_sequence(xp_w[:, K])\n",
    "\n",
    "# Convert results back to assets, consumption and income\n",
    "xp = np.column_stack((xp_w, xp_r[:, 1:]))\n",
    "assets = xp[0, :]                  # Assets\n",
    "\n",
    "up = np.column_stack((up_w, up_r))\n",
    "c = up.flatten() + c_bar           # Consumption\n",
    "\n",
    "time = np.arange(1, K+1)\n",
    "income_w = σ * wp_w[0, 1:K+1] + m1 * time + m2 * time**2  # Income\n",
    "income_r = np.full(T-K, s)\n",
    "income = np.concatenate((income_w, income_r))\n",
    "\n",
    "# Plot results\n",
    "n_rows = 2\n",
    "fig, axes = plt.subplots(n_rows, 1, figsize=(12, 10))\n",
    "\n",
    "plt.subplots_adjust(hspace=0.5)\n",
    "\n",
    "bbox = (0., 1.02, 1., .102)\n",
    "legend_args = {'bbox_to_anchor': bbox, 'loc': 3, 'mode': 'expand'}\n",
    "p_args = {'lw': 2, 'alpha': 0.7}\n",
    "\n",
    "axes[0].plot(range(1, T+1), income, 'g-', label=\"non-financial income\",\n",
    "            **p_args)\n",
    "axes[0].plot(range(T), c, 'k-', label=\"consumption\", **p_args)\n",
    "\n",
    "axes[1].plot(range(T+1), assets, 'b-', label=\"assets\", **p_args)\n",
    "axes[1].plot(range(T+1), np.zeros(T+1), 'k-')\n",
    "\n",
    "for ax in axes:\n",
    "    ax.grid()\n",
    "    ax.set_xlabel('Time')\n",
    "    ax.legend(ncol=2, **legend_args)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0889569",
   "metadata": {},
   "source": [
    "## Exercise 56.3\n",
    "\n",
    "Reproduce the figures from the monopolist application [given above](#lqc-mwac).\n",
    "\n",
    "For parameters, use $ a_0 = 5, a_1 = 0.5, \\sigma = 0.15, \\rho = 0.9,\n",
    "\\beta = 0.95 $ and $ c = 2 $, while $ \\gamma $ varies between 1 and 50\n",
    "(see figures)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61f92342",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 56.3](https://python.quantecon.org/#lqc_ex3)\n",
    "\n",
    "The first task is to find the matrices $ A, B, C, Q, R $ that define\n",
    "the LQ problem.\n",
    "\n",
    "Recall that $ x_t = (\\bar q_t \\;\\, q_t \\;\\, 1)' $, while\n",
    "$ u_t = q_{t+1} - q_t $.\n",
    "\n",
    "Letting $ m_0 := (a_0 - c) / 2a_1 $ and $ m_1 := 1 / 2 a_1 $, we\n",
    "can write $ \\bar q_t = m_0 + m_1 d_t $, and then, with some\n",
    "manipulation\n",
    "\n",
    "$$\n",
    "\\bar q_{t+1} = m_0 (1 - \\rho) + \\rho \\bar q_t + m_1 \\sigma w_{t+1}\n",
    "$$\n",
    "\n",
    "By our definition of $ u_t $, the dynamics of $ q_t $ are\n",
    "$ q_{t+1} = q_t + u_t $.\n",
    "\n",
    "Using these facts you should be able to build the correct\n",
    "$ A, B, C $ matrices (and then check them against those found in the\n",
    "solution code below).\n",
    "\n",
    "Suitable $ R, Q $ matrices can be found by inspecting the objective\n",
    "function, which we repeat here for convenience:\n",
    "\n",
    "$$\n",
    "\\min\n",
    "\\mathbb E \\,\n",
    "\\left\\{\n",
    "    \\sum_{t=0}^{\\infty} \\beta^t\n",
    "    a_1 ( q_t - \\bar q_t)^2 + \\gamma u_t^2\n",
    "\\right\\}\n",
    "$$\n",
    "\n",
    "Our solution code is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0f6f6b6",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Model parameters\n",
    "a0 = 5\n",
    "a1 = 0.5\n",
    "σ = 0.15\n",
    "ρ = 0.9\n",
    "γ = 1\n",
    "β = 0.95\n",
    "c = 2\n",
    "T = 120\n",
    "\n",
    "# Useful constants\n",
    "m0 = (a0-c)/(2 * a1)\n",
    "m1 = 1/(2 * a1)\n",
    "\n",
    "# Formulate LQ problem\n",
    "Q = γ\n",
    "R = [[ a1, -a1,  0],\n",
    "     [-a1,  a1,  0],\n",
    "     [  0,   0,  0]]\n",
    "A = [[ρ, 0, m0 * (1 - ρ)],\n",
    "     [0, 1,            0],\n",
    "     [0, 0,            1]]\n",
    "\n",
    "B = [[0],\n",
    "     [1],\n",
    "     [0]]\n",
    "C = [[m1 * σ],\n",
    "     [     0],\n",
    "     [     0]]\n",
    "\n",
    "lq = LQ(Q, R, A, B, C=C, beta=β)\n",
    "\n",
    "# Simulate state / control paths\n",
    "x0 = (m0, 2, 1)\n",
    "xp, up, wp = lq.compute_sequence(x0, ts_length=150)\n",
    "q_bar = xp[0, :]\n",
    "q = xp[1, :]\n",
    "\n",
    "# Plot simulation results\n",
    "fig, ax = plt.subplots(figsize=(10, 6.5))\n",
    "\n",
    "# Some fancy plotting stuff -- simplify if you prefer\n",
    "bbox = (0., 1.01, 1., .101)\n",
    "legend_args = {'bbox_to_anchor': bbox, 'loc': 3, 'mode': 'expand'}\n",
    "p_args = {'lw': 2, 'alpha': 0.6}\n",
    "\n",
    "time = range(len(q))\n",
    "ax.set(xlabel='Time', xlim=(0, max(time)))\n",
    "ax.plot(time, q_bar, 'k-', lw=2, alpha=0.6, label=r'$\\bar q_t$')\n",
    "ax.plot(time, q, 'b-', lw=2, alpha=0.6, label='$q_t$')\n",
    "ax.legend(ncol=2, **legend_args)\n",
    "s = f'dynamics with $\\gamma = {γ}$'\n",
    "ax.text(max(time) * 0.6, 1 * q_bar.max(), s, fontsize=14)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "date": 1733369121.340165,
  "filename": "lqcontrol.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "LQ Control: Foundations"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}