lqcontrol.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LQ Dynamic Programming Problems\n",
    "\n",
    "\n",
    "<a id='index-0'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [LQ Dynamic Programming Problems](#LQ-Dynamic-Programming-Problems)  \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)  \n",
    "  - [Solutions](#Solutions)  "
   ]
  },
  {
   "cell_type": "markdown",
   "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://lectures.quantecon.org/jl/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://lectures.quantecon.org/jl/lu_tricks.html#) and [this lecture](https://lectures.quantecon.org/jl/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://lectures.quantecon.org/jl/short_path.html#) and [this lecture](https://lectures.quantecon.org/jl/optgrowth.html#))  \n",
    "\n",
    "\n",
    "For additional reading on LQ control, see, for example,\n",
    "\n",
    "- [[LS18]](https://lectures.quantecon.org/jl/zreferences.html#ljungqvist2012), chapter 5  \n",
    "- [[HS08]](https://lectures.quantecon.org/jl/zreferences.html#hansensargent2008), chapter 4  \n",
    "- [[HLL96]](https://lectures.quantecon.org/jl/zreferences.html#hernandezlermalasserre1996), 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)"
   ]
  },
  {
   "cell_type": "markdown",
   "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",
   "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{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",
   "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 [(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",
   "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{2}\n",
    "$$\n",
    "\n",
    "How can we write this new system in the form of equation [(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://lectures.quantecon.org/jl/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}    -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{3}\n",
    "$$\n",
    "\n",
    "then the first row is equivalent to [(2)](#equation-lq-lomwc)\n",
    "\n",
    "Moreover, the model is now linear, and can be written in the form of\n",
    "[(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{4}\n",
    "$$\n",
    "\n",
    "In effect, we’ve bought ourselves linearity by adding another state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Preferences\n",
    "\n",
    "In the LQ model, the aim is to minimize a 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{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 [[HS08]](https://lectures.quantecon.org/jl/zreferences.html#hansensargent2008) for details"
   ]
  },
  {
   "cell_type": "markdown",
   "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 [(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",
   "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",
   "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",
   "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{6}\n",
    "$$\n",
    "\n",
    "subject to the law of motion [(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",
   "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",
   "metadata": {},
   "source": [
    "### Solution\n",
    "\n",
    "To solve the finite horizon LQ problem we can use a dynamic programming\n",
    "strategy based on backwards induction that is conceptually similar to the approach adopted in [this lecture](https://lectures.quantecon.org/jl/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{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 backwards 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, [(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{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://lectures.quantecon.org/jl/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{9}\n",
    "$$\n",
    "\n",
    "Plugging this back into [(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}\n",
    "= 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{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{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}\n",
    "= 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{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{13}\n",
    "$$\n",
    "\n",
    "Recalling [(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{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 [(14)](#equation-lq-oc) and [(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{15}\n",
    "$$\n",
    "\n",
    "for $ t = 0, \\ldots, T-1 $ attains the minimum of [(6)](#equation-lq-object) subject to our constraints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "We will use code from [lqcontrol.jl](https://github.com/QuantEcon/QuantEcon.jl/blob/master/src/lqcontrol.jl)\n",
    "in [QuantEcon.jl](http://quantecon.org/julia_index.html)\n",
    "to solve finite and infinite horizon linear quadratic control problems\n",
    "\n",
    "In the module, the various updating, simulation and fixed point methods  act on a type  called `LQ`, which includes\n",
    "\n",
    "- Instance data:  \n",
    "  \n",
    "  - The required parameters $ Q, R, A, B $ and optional parameters C, β, T, R_f, N specifying a given LQ model  \n",
    "    \n",
    "    - set $ T $ and $ R_f $ to `None` in the infinite horizon case  \n",
    "    - set `C = None` (or zero) in the deterministic case  \n",
    "    \n",
    "  - the value function and policy data  \n",
    "    \n",
    "    - $ d_t, P_t, F_t $ in the finite horizon case  \n",
    "    - $ d, P, F $ in the infinite horizon case  \n",
    "    \n",
    "  \n",
    "- Methods:  \n",
    "  \n",
    "  - `update_values` — shifts $ d_t, P_t, F_t $ to their $ t-1 $ values via [(12)](#equation-lq-pr), [(13)](#equation-lq-dd) and [(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",
    "\n",
    "<a id='lq-mfpa'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "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, [[Fri56]](https://lectures.quantecon.org/jl/zreferences.html#friedman1956) or [[MB54]](https://lectures.quantecon.org/jl/zreferences.html#modiglianibrumberg1954))\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{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 [(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 [(4)](#equation-lq-lowmc2)\n",
    "\n",
    "To match with this state and control, the objective function [(16)](#equation-lq-pio) can\n",
    "be written in the form of [(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 [(12)](#equation-lq-pr) and [(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 [(15)](#equation-lq-xud)\n",
    "\n",
    "The following figure was computed using $ r = 0.05, \\beta = 1 / (1 + r), \\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": "markdown",
   "metadata": {},
   "source": [
    "### Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "using InstantiateFromURL\n",
    "\n",
    "# activate the QuantEcon environment\n",
    "activate_github(\"QuantEcon/QuantEconLecturePackages\", tag = \"v0.9.5\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "using LinearAlgebra, Statistics, Compat\n",
    "using Plots, Plots.PlotMeasures, QuantEcon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"600\" viewBox=\"0 0 2400 2400\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip7100\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2000\" height=\"2000\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7101\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"2400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip7101)\" points=\"\n",
       "0,2400 2400,2400 2400,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7102\">\n",
       "    <rect x=\"480\" y=\"240\" width=\"1681\" height=\"1681\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip7101)\" points=\"\n",
       "189.504,741.27 2321.26,741.27 2321.26,47.2441 189.504,47.2441 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7103\">\n",
       "    <rect x=\"189\" y=\"47\" width=\"2133\" height=\"695\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip7103)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  205.146,741.27 205.146,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7103)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  652.055,741.27 652.055,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7103)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1098.96,741.27 1098.96,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7103)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1545.87,741.27 1545.87,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7103)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1992.78,741.27 1992.78,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7103)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,645.515 2321.26,645.515 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7103)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,498.809 2321.26,498.809 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7103)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,352.102 2321.26,352.102 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7103)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,205.396 2321.26,205.396 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7103)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,58.6897 2321.26,58.6897 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,741.27 2321.26,741.27 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,741.27 189.504,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  205.146,741.27 205.146,730.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  652.055,741.27 652.055,730.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1098.96,741.27 1098.96,730.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1545.87,741.27 1545.87,730.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1992.78,741.27 1992.78,730.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,645.515 221.48,645.515 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,498.809 221.48,498.809 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,352.102 221.48,352.102 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,205.396 221.48,205.396 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,58.6897 221.48,58.6897 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 205.146, 795.27)\" x=\"205.146\" y=\"795.27\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 652.055, 795.27)\" x=\"652.055\" y=\"795.27\">10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1098.96, 795.27)\" x=\"1098.96\" y=\"795.27\">20</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1545.87, 795.27)\" x=\"1545.87\" y=\"795.27\">30</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1992.78, 795.27)\" x=\"1992.78\" y=\"795.27\">40</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 663.015)\" x=\"165.504\" y=\"663.015\">-1.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 516.309)\" x=\"165.504\" y=\"516.309\">-0.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 369.602)\" x=\"165.504\" y=\"369.602\">0.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 222.896)\" x=\"165.504\" y=\"222.896\">0.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 76.1897)\" x=\"165.504\" y=\"76.1897\">1.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(0, 1255.38, 1190.4)\" x=\"1255.38\" y=\"1190.4\">Time</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip7103)\" style=\"stroke:#0000ff; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  249.837,352.102 294.528,416.97 339.219,464.559 383.91,501.99 428.6,462.083 473.291,508.046 517.982,477.088 562.673,584.528 607.364,628.435 652.055,721.628 \n",
       "  696.746,717.54 741.437,649.993 786.128,588.111 830.818,575.503 875.509,489.383 920.2,407.376 964.891,375.758 1009.58,520.999 1054.27,378.582 1098.96,317.602 \n",
       "  1143.65,257.817 1188.35,261.252 1233.04,312.492 1277.73,303.492 1322.42,215.308 1367.11,204.819 1411.8,203.542 1456.49,66.8863 1501.18,81.1216 1545.87,154.719 \n",
       "  1590.56,219.519 1635.25,181.582 1679.95,110.97 1724.64,204.381 1769.33,272.666 1814.02,321.289 1858.71,426.833 1903.4,469.574 1948.09,524.741 1992.78,516.185 \n",
       "  2037.47,463.204 2082.16,497.087 2126.85,468.89 2171.55,434.167 2216.24,494.772 2260.93,208.18 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7103)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  249.837,352.102 294.528,352.102 339.219,352.102 383.91,352.102 428.6,352.102 473.291,352.102 517.982,352.102 562.673,352.102 607.364,352.102 652.055,352.102 \n",
       "  696.746,352.102 741.437,352.102 786.128,352.102 830.818,352.102 875.509,352.102 920.2,352.102 964.891,352.102 1009.58,352.102 1054.27,352.102 1098.96,352.102 \n",
       "  1143.65,352.102 1188.35,352.102 1233.04,352.102 1277.73,352.102 1322.42,352.102 1367.11,352.102 1411.8,352.102 1456.49,352.102 1501.18,352.102 1545.87,352.102 \n",
       "  1590.56,352.102 1635.25,352.102 1679.95,352.102 1724.64,352.102 1769.33,352.102 1814.02,352.102 1858.71,352.102 1903.4,352.102 1948.09,352.102 1992.78,352.102 \n",
       "  2037.47,352.102 2082.16,352.102 2126.85,352.102 2171.55,352.102 2216.24,352.102 2260.93,352.102 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7103)\" style=\"stroke:#ff0000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  249.837,368.319 294.528,380.324 339.219,389.878 383.91,380.179 428.6,391.884 473.291,384.468 517.982,411.601 562.673,423.118 607.364,447.096 652.055,447.038 \n",
       "  696.746,431.162 741.437,416.567 786.128,414.152 830.818,393.364 875.509,373.347 920.2,365.651 964.891,402.056 1009.58,367.175 1054.27,352.051 1098.96,336.936 \n",
       "  1143.65,337.301 1188.35,349.6 1233.04,347.111 1277.73,324.75 1322.42,321.17 1367.11,319.737 1411.8,284.357 1456.49,285.381 1501.18,301.159 1545.87,315.273 \n",
       "  1590.56,304.253 1635.25,284.425 1679.95,304.374 1724.64,319.125 1769.33,329.883 1814.02,355.656 1858.71,368.036 1903.4,384.903 1948.09,388.065 1992.78,380.851 \n",
       "  2037.47,394.348 2082.16,395.708 2126.85,396.289 2171.55,421.448 2216.24,385.467 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip7101)\" points=\"\n",
       "1311.03,312.204 2249.26,312.204 2249.26,130.764 1311.03,130.764 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1311.03,312.204 2249.26,312.204 2249.26,130.764 1311.03,130.764 1311.03,312.204 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#0000ff; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1335.03,191.244 1479.03,191.244 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1503.03, 208.744)\" x=\"1503.03\" y=\"208.744\">assets</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#ff0000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1335.03,251.724 1479.03,251.724 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1503.03, 269.224)\" x=\"1503.03\" y=\"269.224\">cumulative unanticipated income</text>\n",
       "</g>\n",
       "<polygon clip-path=\"url(#clip7101)\" points=\"\n",
       "189.504,1941.27 2321.26,1941.27 2321.26,1247.24 189.504,1247.24 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7104\">\n",
       "    <rect x=\"189\" y=\"1247\" width=\"2133\" height=\"695\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip7104)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  204.13,1941.27 204.13,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7104)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  661.196,1941.27 661.196,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7104)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1118.26,1941.27 1118.26,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7104)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1575.33,1941.27 1575.33,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7104)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  2032.39,1941.27 2032.39,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7104)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,1910.72 2321.26,1910.72 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7104)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,1807.31 2321.26,1807.31 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7104)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,1703.91 2321.26,1703.91 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7104)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,1600.5 2321.26,1600.5 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7104)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,1497.1 2321.26,1497.1 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7104)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,1393.69 2321.26,1393.69 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7104)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,1290.29 2321.26,1290.29 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1941.27 2321.26,1941.27 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1941.27 189.504,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  204.13,1941.27 204.13,1930.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  661.196,1941.27 661.196,1930.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1118.26,1941.27 1118.26,1930.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1575.33,1941.27 1575.33,1930.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2032.39,1941.27 2032.39,1930.86 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1910.72 221.48,1910.72 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1807.31 221.48,1807.31 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1703.91 221.48,1703.91 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1600.5 221.48,1600.5 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1497.1 221.48,1497.1 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1393.69 221.48,1393.69 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1290.29 221.48,1290.29 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 204.13, 1995.27)\" x=\"204.13\" y=\"1995.27\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 661.196, 1995.27)\" x=\"661.196\" y=\"1995.27\">10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1118.26, 1995.27)\" x=\"1118.26\" y=\"1995.27\">20</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1575.33, 1995.27)\" x=\"1575.33\" y=\"1995.27\">30</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2032.39, 1995.27)\" x=\"2032.39\" y=\"1995.27\">40</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 1928.22)\" x=\"165.504\" y=\"1928.22\">0.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 1824.81)\" x=\"165.504\" y=\"1824.81\">0.6</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 1721.41)\" x=\"165.504\" y=\"1721.41\">0.7</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 1618)\" x=\"165.504\" y=\"1618\">0.8</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 1514.6)\" x=\"165.504\" y=\"1514.6\">0.9</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 1411.19)\" x=\"165.504\" y=\"1411.19\">1.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 1307.79)\" x=\"165.504\" y=\"1307.79\">1.1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(0, 1255.38, 2390.4)\" x=\"1255.38\" y=\"2390.4\">Time</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip7104)\" style=\"stroke:#008000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  249.837,1393.69 295.543,1406.63 341.25,1416.28 386.957,1424.01 432.663,1416.1 478.37,1425.72 524.076,1419.58 569.783,1442.25 615.49,1451.97 661.196,1472.39 \n",
       "  706.903,1472.34 752.609,1458.52 798.316,1445.67 844.023,1443.51 889.729,1424.72 935.436,1406.36 981.142,1399.2 1026.85,1433.64 1072.56,1400.06 1118.26,1385.23 \n",
       "  1163.97,1370.11 1209.68,1370.49 1255.38,1383.34 1301.09,1380.68 1346.8,1356.09 1392.5,1352.04 1438.21,1350.37 1483.91,1307.7 1529.62,1308.98 1575.33,1329.51 \n",
       "  1621.03,1348.67 1666.74,1332.98 1712.45,1303.22 1758.15,1334.95 1803.86,1359.99 1849.57,1379.63 1895.27,1430.74 1940.98,1457.74 1986.69,1498.84 2032.39,1507.62 \n",
       "  2078.1,1484.13 2123.81,1537.79 2169.51,1544.83 2215.22,1549.23 2260.93,1921.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7104)\" style=\"stroke:#ffa500; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  249.837,1450.84 295.543,1436 341.25,1427.36 386.957,1359.51 432.663,1434.94 478.37,1367.55 524.076,1489.32 569.783,1434.28 615.49,1478.19 661.196,1393.49 \n",
       "  706.903,1337.74 752.609,1342.25 798.316,1385.18 844.023,1320.43 889.729,1323.15 935.436,1366.57 981.142,1521.99 1026.85,1270.76 1072.56,1340.39 1118.26,1340.42 \n",
       "  1163.97,1394.98 1209.68,1437.04 1255.38,1384.92 1301.09,1314.88 1346.8,1381.08 1392.5,1388.64 1438.21,1269 1483.91,1397.3 1529.62,1449.29 1575.33,1443.43 \n",
       "  1621.03,1354.85 1666.74,1323.81 1712.45,1464 1758.15,1445.68 1803.86,1431.6 1849.57,1484.52 1895.27,1437.32 1940.98,1453.13 1986.69,1404.83 2032.39,1368.27 \n",
       "  2078.1,1441.26 2123.81,1398.48 2169.51,1395.74 2215.22,1482.36 2260.93,1266.89 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip7101)\" points=\"\n",
       "1555.79,1512.2 2249.26,1512.2 2249.26,1330.76 1555.79,1330.76 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1555.79,1512.2 2249.26,1512.2 2249.26,1330.76 1555.79,1330.76 1555.79,1512.2 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#008000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1579.79,1391.24 1723.79,1391.24 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1747.79, 1408.74)\" x=\"1747.79\" y=\"1408.74\">consumption</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip7101)\" style=\"stroke:#ffa500; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1579.79,1451.72 1723.79,1451.72 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7101)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1747.79, 1469.22)\" x=\"1747.79\" y=\"1469.22\">non-financial income</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# model parameters\n",
    "r = 0.05\n",
    "β = 1 / (1 + r)\n",
    "T = 45\n",
    "c_bar = 2.0\n",
    "σ = 0.25\n",
    "μ = 1.0\n",
    "q = 1e6\n",
    "\n",
    "# formulate as an LQ problem\n",
    "Q = 1.0\n",
    "R = zeros(2, 2)\n",
    "Rf = zeros(2, 2); Rf[1, 1] = q\n",
    "A = [1 + r -c_bar + μ; 0  1]\n",
    "B = [-1.0, 0]\n",
    "C = [σ, 0]\n",
    "\n",
    "# compute solutions and simulate\n",
    "lq = QuantEcon.LQ(Q, R, A, B, C; bet = β, capT = T, rf = Rf)\n",
    "x0 = [0.0, 1]\n",
    "xp, up, wp = compute_sequence(lq, x0)\n",
    "\n",
    "# convert back to assets, consumption and income\n",
    "assets = vec(xp[1, :]) # a_t\n",
    "c = vec(up .+ c_bar) # c_t\n",
    "income = vec(σ * wp[1, 2:end] .+ μ) # y_t\n",
    "\n",
    "# plot results\n",
    "p = plot([assets, c, zeros(T + 1), income, cumsum(income .- μ)],\n",
    "         lab = [\"assets\" \"consumption\" \"\" \"non-financial income\" \"cumulative unanticipated income\"],\n",
    "         color = [:blue :green :black :orange :red],\n",
    "         xaxis = \"Time\", layout = (2, 1),\n",
    "         bottom_margin = 20mm, size = (600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "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": 4,
   "metadata": {
    "hide-output": false,
    "html-class": "collapse"
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"600\" viewBox=\"0 0 2400 2400\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip7300\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2000\" height=\"2000\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7301\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"2400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip7301)\" points=\"\n",
       "0,2400 2400,2400 2400,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7302\">\n",
       "    <rect x=\"480\" y=\"240\" width=\"1681\" height=\"1681\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip7301)\" points=\"\n",
       "161.394,724.503 2321.26,724.503 2321.26,47.2441 161.394,47.2441 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7303\">\n",
       "    <rect x=\"161\" y=\"47\" width=\"2161\" height=\"678\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip7303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  177.242,724.503 177.242,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  630.044,724.503 630.044,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1082.85,724.503 1082.85,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1535.65,724.503 1535.65,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1988.45,724.503 1988.45,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  161.394,635.625 2321.26,635.625 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  161.394,478.515 2321.26,478.515 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  161.394,321.404 2321.26,321.404 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7303)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  161.394,164.294 2321.26,164.294 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,724.503 2321.26,724.503 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,724.503 161.394,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  177.242,724.503 177.242,714.345 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  630.044,724.503 630.044,714.345 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1082.85,724.503 1082.85,714.345 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1535.65,724.503 1535.65,714.345 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1988.45,724.503 1988.45,714.345 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,635.625 193.792,635.625 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,478.515 193.792,478.515 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,321.404 193.792,321.404 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,164.294 193.792,164.294 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 177.242, 778.503)\" x=\"177.242\" y=\"778.503\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 630.044, 778.503)\" x=\"630.044\" y=\"778.503\">10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1082.85, 778.503)\" x=\"1082.85\" y=\"778.503\">20</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1535.65, 778.503)\" x=\"1535.65\" y=\"778.503\">30</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1988.45, 778.503)\" x=\"1988.45\" y=\"778.503\">40</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 137.394, 653.125)\" x=\"137.394\" y=\"653.125\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 137.394, 496.015)\" x=\"137.394\" y=\"496.015\">1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 137.394, 338.904)\" x=\"137.394\" y=\"338.904\">2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 137.394, 181.794)\" x=\"137.394\" y=\"181.794\">3</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(0, 1241.33, 1190.4)\" x=\"1241.33\" y=\"1190.4\">Time</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip7303)\" style=\"stroke:#0000ff; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  222.522,635.625 267.803,640.355 313.083,582.245 358.363,594.092 403.643,538.3 448.923,474.321 494.204,473.288 539.484,388.622 584.764,284.655 630.044,332.872 \n",
       "  675.324,287.842 720.605,300.87 765.885,325.265 811.165,315.467 856.445,255.403 901.725,168.632 947.006,182.97 992.286,154.796 1037.57,108.469 1082.85,81.5742 \n",
       "  1128.13,88.7287 1173.41,66.4118 1218.69,112.9 1263.97,225.848 1309.25,248.414 1354.53,334.495 1399.81,443.667 1445.09,477.071 1490.37,543.658 1535.65,483.414 \n",
       "  1580.93,522.925 1626.21,598.023 1671.49,563.486 1716.77,613.762 1762.05,653.357 1807.33,688.931 1852.61,671.437 1897.89,677.596 1943.17,636.944 1988.45,644.178 \n",
       "  2033.73,622.702 2079.01,621.562 2124.29,619.469 2169.57,605.16 2214.85,599.396 2260.13,628.956 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7303)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  222.522,635.625 267.803,635.625 313.083,635.625 358.363,635.625 403.643,635.625 448.923,635.625 494.204,635.625 539.484,635.625 584.764,635.625 630.044,635.625 \n",
       "  675.324,635.625 720.605,635.625 765.885,635.625 811.165,635.625 856.445,635.625 901.725,635.625 947.006,635.625 992.286,635.625 1037.57,635.625 1082.85,635.625 \n",
       "  1128.13,635.625 1173.41,635.625 1218.69,635.625 1263.97,635.625 1309.25,635.625 1354.53,635.625 1399.81,635.625 1445.09,635.625 1490.37,635.625 1535.65,635.625 \n",
       "  1580.93,635.625 1626.21,635.625 1671.49,635.625 1716.77,635.625 1762.05,635.625 1807.33,635.625 1852.61,635.625 1897.89,635.625 1943.17,635.625 1988.45,635.625 \n",
       "  2033.73,635.625 2079.01,635.625 2124.29,635.625 2169.57,635.625 2214.85,635.625 2260.13,635.625 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7303)\" style=\"stroke:#ff0000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  222.522,641.387 267.803,631.396 313.083,638.638 358.363,628.944 403.643,616.937 448.923,620.356 494.204,602.778 539.484,579.957 584.764,594.662 630.044,586.134 \n",
       "  675.324,591.808 720.605,600.248 765.885,600.116 811.165,587.232 856.445,567.222 901.725,571.87 947.006,565.771 992.286,554.767 1037.57,548.118 1082.85,549.559 \n",
       "  1128.13,543.398 1173.41,553.985 1218.69,581.2 1263.97,586.368 1309.25,607.294 1354.53,634.435 1399.81,643.372 1445.09,660.699 1490.37,646.806 1535.65,656.944 \n",
       "  1580.93,676.211 1626.21,668.862 1671.49,682.014 1716.77,693.101 1762.05,703.752 1807.33,701.788 1852.61,705.336 1897.89,697.336 1943.17,700.069 1988.45,695.763 \n",
       "  2033.73,695.48 2079.01,694.581 2124.29,690.014 2169.57,684.994 2214.85,683.327 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip7301)\" points=\"\n",
       "1311.03,312.204 2249.26,312.204 2249.26,130.764 1311.03,130.764 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1311.03,312.204 2249.26,312.204 2249.26,130.764 1311.03,130.764 1311.03,312.204 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#0000ff; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1335.03,191.244 1479.03,191.244 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1503.03, 208.744)\" x=\"1503.03\" y=\"208.744\">assets</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#ff0000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1335.03,251.724 1479.03,251.724 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1503.03, 269.224)\" x=\"1503.03\" y=\"269.224\">cumulative unanticipated income</text>\n",
       "</g>\n",
       "<polygon clip-path=\"url(#clip7301)\" points=\"\n",
       "161.394,1924.5 2321.26,1924.5 2321.26,1247.24 161.394,1247.24 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7304\">\n",
       "    <rect x=\"161\" y=\"1247\" width=\"2161\" height=\"678\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip7304)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  176.213,1924.5 176.213,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7304)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  639.306,1924.5 639.306,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7304)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1102.4,1924.5 1102.4,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7304)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1565.49,1924.5 1565.49,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7304)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  2028.59,1924.5 2028.59,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7304)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  161.394,1792.7 2321.26,1792.7 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7304)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  161.394,1638.87 2321.26,1638.87 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7304)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  161.394,1485.04 2321.26,1485.04 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7304)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  161.394,1331.21 2321.26,1331.21 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,1924.5 2321.26,1924.5 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,1924.5 161.394,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  176.213,1924.5 176.213,1914.34 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  639.306,1924.5 639.306,1914.34 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1102.4,1924.5 1102.4,1914.34 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1565.49,1924.5 1565.49,1914.34 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2028.59,1924.5 2028.59,1914.34 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,1792.7 193.792,1792.7 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,1638.87 193.792,1638.87 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,1485.04 193.792,1485.04 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,1331.21 193.792,1331.21 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 176.213, 1978.5)\" x=\"176.213\" y=\"1978.5\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 639.306, 1978.5)\" x=\"639.306\" y=\"1978.5\">10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1102.4, 1978.5)\" x=\"1102.4\" y=\"1978.5\">20</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1565.49, 1978.5)\" x=\"1565.49\" y=\"1978.5\">30</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2028.59, 1978.5)\" x=\"2028.59\" y=\"1978.5\">40</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 137.394, 1810.2)\" x=\"137.394\" y=\"1810.2\">0.9</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 137.394, 1656.37)\" x=\"137.394\" y=\"1656.37\">1.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 137.394, 1502.54)\" x=\"137.394\" y=\"1502.54\">1.1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 137.394, 1348.71)\" x=\"137.394\" y=\"1348.71\">1.2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(0, 1241.33, 2390.4)\" x=\"1241.33\" y=\"2390.4\">Time</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip7304)\" style=\"stroke:#008000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  222.522,1818.22 268.832,1818.84 315.141,1780.37 361.45,1785.14 407.76,1747.39 454.069,1703.94 500.378,1699.93 546.687,1642.13 592.997,1570.89 639.306,1597.55 \n",
       "  685.615,1563.25 731.925,1566.7 778.234,1577.69 824.543,1565.61 870.853,1518.67 917.162,1451.9 963.471,1454.24 1009.78,1426.11 1056.09,1383.74 1102.4,1354 \n",
       "  1148.71,1348.36 1195.02,1319.67 1241.33,1342.95 1287.64,1419.68 1333.95,1426.05 1380.25,1485.5 1426.56,1567.84 1472.87,1587.88 1519.18,1639.59 1565.49,1574.55 \n",
       "  1611.8,1602.99 1658.11,1670.8 1704.42,1626.39 1750.73,1674.7 1797.04,1716.37 1843.35,1759.35 1889.66,1735.05 1935.97,1744.13 1982.28,1675.72 2028.59,1684.71 \n",
       "  2074.89,1632.6 2121.2,1617.27 2167.51,1592.2 2213.82,1483.8 2260.13,1266.41 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7304)\" style=\"stroke:#ffa500; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  222.522,1695.28 268.832,1541.04 315.141,1709.78 361.45,1543.95 407.76,1521.31 454.069,1672.35 500.378,1466.76 546.687,1415.43 592.997,1782.85 639.306,1555.37 \n",
       "  685.615,1694.42 731.925,1721.51 778.234,1637.58 824.543,1512.72 870.853,1442.96 917.162,1684.38 963.471,1579.15 1009.78,1531.13 1056.09,1573.77 1102.4,1652.97 \n",
       "  1148.71,1578.55 1195.02,1742.53 1241.33,1905.34 1287.64,1689.46 1333.95,1843.76 1380.25,1904.61 1426.56,1726.37 1472.87,1808.52 1519.18,1502.84 1565.49,1738.13 \n",
       "  1611.8,1827.51 1658.11,1566.92 1704.42,1767.64 1750.73,1747.42 1797.04,1743.15 1843.35,1619.64 1889.66,1673.61 1935.97,1560.54 1982.28,1665.63 2028.59,1596.71 \n",
       "  2074.89,1636.09 2121.2,1630.07 2167.51,1594.16 2213.82,1589.72 2260.13,1622.55 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip7301)\" points=\"\n",
       "1555.79,1512.2 2249.26,1512.2 2249.26,1330.76 1555.79,1330.76 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1555.79,1512.2 2249.26,1512.2 2249.26,1330.76 1555.79,1330.76 1555.79,1512.2 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#008000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1579.79,1391.24 1723.79,1391.24 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1747.79, 1408.74)\" x=\"1747.79\" y=\"1408.74\">consumption</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip7301)\" style=\"stroke:#ffa500; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1579.79,1451.72 1723.79,1451.72 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7301)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1747.79, 1469.22)\" x=\"1747.79\" y=\"1469.22\">non-financial income</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute solutions and simulate\n",
    "lq = QuantEcon.LQ(Q, R, A, B, C; bet = 0.96, capT = T, rf = Rf)\n",
    "x0 = [0.0, 1]\n",
    "xp, up, wp = compute_sequence(lq, x0)\n",
    "\n",
    "# convert back to assets, consumption and income\n",
    "assets = vec(xp[1, :]) # a_t\n",
    "c = vec(up .+ c_bar) # c_t\n",
    "income = vec(σ * wp[1, 2:end] .+ μ) # y_t\n",
    "\n",
    "# plot results\n",
    "p = plot([assets, c, zeros(T + 1), income, cumsum(income .- μ)],\n",
    "         lab = [\"assets\" \"consumption\" \"\" \"non-financial income\" \"cumulative unanticipated income\"],\n",
    "         color = [:blue :green :black :orange :red],\n",
    "         xaxis = \"Time\", layout = (2, 1),\n",
    "         bottom_margin = 20mm, size = (600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have a slowly rising consumption stream and a hump-shaped build\n",
    "up 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",
   "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",
   "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 [[HS13]](https://lectures.quantecon.org/jl/zreferences.html#hansensargent2013), section 2.4\n",
    "\n",
    "\n",
    "<a id='lq-cpt'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "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{17}\n",
    "$$\n",
    "\n",
    "Our results extend to this case in a straightforward way\n",
    "\n",
    "The sequence $ \\{P_t\\} $ from [(12)](#equation-lq-pr) becomes\n",
    "\n",
    "\n",
    "<a id='equation-lq-pr-cp'></a>\n",
    "$$\n",
    "P_{t-1}\n",
    "= R - (\\beta B' P_t A + 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{18}\n",
    "$$\n",
    "\n",
    "The policies in [(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{19}\n",
    "$$\n",
    "\n",
    "The sequence $ \\{d_t\\} $ is unchanged from [(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",
   "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{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)' (Q + \\beta B' P B)^{-1} (\\beta B' P A + N) +\n",
    "\\beta A' P A \\tag{21}\n",
    "$$\n",
    "\n",
    "Equation [(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 [(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{22}\n",
    "$$\n",
    "\n",
    "The sequence $ \\{d_t\\} $ from [(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{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",
   "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 [(22)](#equation-lq-oc-ih) or [(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",
   "metadata": {},
   "source": [
    "## Further Applications\n",
    "\n",
    "\n",
    "<a id='lq-nsi'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "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 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{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{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{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 [(25)](#equation-lq-hib) as desired,\n",
    "and that the other state variables also update appropriately\n",
    "\n",
    "To implement preference specification [(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{27}\n",
    "$$\n",
    "\n",
    "The next figure shows a simulation of consumption and assets computed using\n",
    "the `compute_sequence` method of `lqcontrol.jl` with initial assets set to zero\n",
    "\n",
    "\n",
    "<a id='solution-lqc-ex1-fig'></a>\n",
    "<img src=\"https://s3-ap-southeast-2.amazonaws.com/lectures.quantecon.org/jl/_static/figures/solution_lqc_ex1.png\" style=\"\">\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 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",
   "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{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 [(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 [(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 backwards 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 backwards 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",
    "<img src=\"https://s3-ap-southeast-2.amazonaws.com/lectures.quantecon.org/jl/_static/figures/solution_lqc_ex2.png\" style=\"\">\n",
    "\n",
    "  \n",
    "The full set of parameters used in the simulation is discussed in [Exercise 2](#lqc-ex2), 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",
   "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{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",
    "<img src=\"https://s3-ap-southeast-2.amazonaws.com/lectures.quantecon.org/jl/_static/figures/solution_lqc_ex3_g1.png\" style=\"\">\n",
    "\n",
    "  \n",
    "<img src=\"https://s3-ap-southeast-2.amazonaws.com/lectures.quantecon.org/jl/_static/figures/solution_lqc_ex3_g10.png\" style=\"\">\n",
    "\n",
    "  \n",
    "<img src=\"https://s3-ap-southeast-2.amazonaws.com/lectures.quantecon.org/jl/_static/figures/solution_lqc_ex3_g50.png\" style=\"\">\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 [(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 [(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 [(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{30}\n",
    "$$\n",
    "\n",
    "It’s now relatively straightforward to find $ R $ and $ Q $ such that\n",
    "[(30)](#equation-lq-object-mp2) can be written as [(20)](#equation-lq-object-ih)\n",
    "\n",
    "Furthermore, the matrices $ A, B $ and $ C $ from [(1)](#equation-lq-lom)\n",
    "can be found by writing down the dynamics of each element of the state\n",
    "\n",
    "[Exercise 3](#lqc-ex3) asks you to complete this process, and reproduce the preceding figures"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "\n",
    "<a id='lqc-ex1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 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 $\n",
    "\n",
    "\n",
    "<a id='lqc-ex2'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 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",
    "Some hints are as follows:\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 [(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 [(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\n",
    "\n",
    "\n",
    "<a id='lqc-ex3'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 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",
   "metadata": {},
   "source": [
    "## Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\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": 5,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"600\" viewBox=\"0 0 2400 2400\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip7500\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2000\" height=\"2000\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7501\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"2400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip7501)\" points=\"\n",
       "0,2400 2400,2400 2400,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7502\">\n",
       "    <rect x=\"480\" y=\"240\" width=\"1681\" height=\"1681\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip7501)\" points=\"\n",
       "189.504,724.503 2321.26,724.503 2321.26,47.2441 189.504,47.2441 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7503\">\n",
       "    <rect x=\"189\" y=\"47\" width=\"2133\" height=\"678\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip7503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  208.794,724.503 208.794,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  619.221,724.503 619.221,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1029.65,724.503 1029.65,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1440.07,724.503 1440.07,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1850.5,724.503 1850.5,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  2260.93,724.503 2260.93,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,704.092 2321.26,704.092 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,547.604 2321.26,547.604 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,391.116 2321.26,391.116 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,234.628 2321.26,234.628 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7503)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,78.14 2321.26,78.14 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,724.503 2321.26,724.503 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,724.503 189.504,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  208.794,724.503 208.794,714.345 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  619.221,724.503 619.221,714.345 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1029.65,724.503 1029.65,714.345 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1440.07,724.503 1440.07,714.345 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1850.5,724.503 1850.5,714.345 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2260.93,724.503 2260.93,714.345 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,704.092 221.48,704.092 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,547.604 221.48,547.604 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,391.116 221.48,391.116 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,234.628 221.48,234.628 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,78.14 221.48,78.14 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 208.794, 778.503)\" x=\"208.794\" y=\"778.503\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 619.221, 778.503)\" x=\"619.221\" y=\"778.503\">10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1029.65, 778.503)\" x=\"1029.65\" y=\"778.503\">20</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1440.07, 778.503)\" x=\"1440.07\" y=\"778.503\">30</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1850.5, 778.503)\" x=\"1850.5\" y=\"778.503\">40</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2260.93, 778.503)\" x=\"2260.93\" y=\"778.503\">50</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 721.592)\" x=\"165.504\" y=\"721.592\">0.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 565.104)\" x=\"165.504\" y=\"565.104\">0.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 408.616)\" x=\"165.504\" y=\"408.616\">1.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 252.128)\" x=\"165.504\" y=\"252.128\">1.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 95.64)\" x=\"165.504\" y=\"95.64\">2.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(0, 1255.38, 1190.4)\" x=\"1255.38\" y=\"1190.4\">Time</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip7503)\" style=\"stroke:#ffa500; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  249.837,653.886 290.879,608.077 331.922,554.396 372.965,523.042 414.007,482.877 455.05,439.299 496.093,405.9 537.135,374.696 578.178,340.669 619.221,303.324 \n",
       "  660.263,276.326 701.306,240.663 742.349,214.946 783.391,203.797 824.434,164.354 865.477,149.106 906.519,152.671 947.562,142.471 988.605,111.831 1029.65,95.1466 \n",
       "  1070.69,88.0118 1111.73,75.5599 1152.78,92.3954 1193.82,66.4118 1234.86,88.8909 1275.9,67.9633 1316.95,71.5763 1357.99,91.9551 1399.03,93.7342 1440.07,95.1153 \n",
       "  1481.12,115.379 1522.16,125.647 1563.2,134.523 1604.24,152.92 1645.29,176.397 1686.33,194.345 1727.37,231.709 1768.42,245.287 1809.46,264.591 1850.5,300.006 \n",
       "  1891.54,324.497 1932.59,363.768 1973.63,393.421 2014.67,442.259 2055.71,485.946 2096.76,524.298 2137.8,562.467 2178.84,608.985 2219.88,646.466 2260.93,705.336 \n",
       "  \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7503)\" style=\"stroke:#008000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  249.837,332.458 290.879,332.042 331.922,332.09 372.965,328.946 414.007,330.15 455.05,331.698 496.093,331.551 537.135,332.792 578.178,335.514 619.221,337.881 \n",
       "  660.263,337.819 701.306,338.558 742.349,335.897 783.391,332.939 824.434,334.741 865.477,329.066 906.519,324.884 947.562,329.231 988.605,335.667 1029.65,334.656 \n",
       "  1070.69,331.173 1111.73,328.464 1152.78,323.276 1193.82,327.942 1234.86,322.039 1275.9,327.124 1316.95,321.723 1357.99,316.499 1399.03,318.931 1440.07,318.707 \n",
       "  1481.12,314.394 1522.16,315.047 1563.2,314.153 1604.24,309.591 1645.29,305.689 1686.33,304.472 1727.37,301.118 1768.42,307.753 1809.46,306.166 1850.5,298.262 \n",
       "  1891.54,295.26 1932.59,285.85 1973.63,281.917 2014.67,271.303 2055.71,274.677 2096.76,285.758 2137.8,294.192 2178.84,293.188 2219.88,296.911 2260.93,237.053 \n",
       "  \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip7501)\" points=\"\n",
       "1555.79,312.204 2249.26,312.204 2249.26,130.764 1555.79,130.764 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1555.79,312.204 2249.26,312.204 2249.26,130.764 1555.79,130.764 1555.79,312.204 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#ffa500; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1579.79,191.244 1723.79,191.244 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1747.79, 208.744)\" x=\"1747.79\" y=\"208.744\">non-financial income</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#008000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1579.79,251.724 1723.79,251.724 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1747.79, 269.224)\" x=\"1747.79\" y=\"269.224\">consumption</text>\n",
       "</g>\n",
       "<polygon clip-path=\"url(#clip7501)\" points=\"\n",
       "189.504,1924.5 2321.26,1924.5 2321.26,1247.24 189.504,1247.24 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7504\">\n",
       "    <rect x=\"189\" y=\"1247\" width=\"2133\" height=\"678\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip7504)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  209.615,1924.5 209.615,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7504)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  611.833,1924.5 611.833,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7504)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1014.05,1924.5 1014.05,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7504)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1416.27,1924.5 1416.27,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7504)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1818.49,1924.5 1818.49,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7504)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  2220.71,1924.5 2220.71,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7504)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,1855.36 2321.26,1855.36 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7504)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,1732.84 2321.26,1732.84 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7504)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,1610.32 2321.26,1610.32 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7504)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,1487.8 2321.26,1487.8 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7504)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  189.504,1365.28 2321.26,1365.28 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1924.5 2321.26,1924.5 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1924.5 189.504,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  209.615,1924.5 209.615,1914.34 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  611.833,1924.5 611.833,1914.34 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1014.05,1924.5 1014.05,1914.34 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1416.27,1924.5 1416.27,1914.34 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1818.49,1924.5 1818.49,1914.34 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2220.71,1924.5 2220.71,1914.34 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1855.36 221.48,1855.36 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1732.84 221.48,1732.84 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1610.32 221.48,1610.32 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1487.8 221.48,1487.8 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  189.504,1365.28 221.48,1365.28 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 209.615, 1978.5)\" x=\"209.615\" y=\"1978.5\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 611.833, 1978.5)\" x=\"611.833\" y=\"1978.5\">10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1014.05, 1978.5)\" x=\"1014.05\" y=\"1978.5\">20</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1416.27, 1978.5)\" x=\"1416.27\" y=\"1978.5\">30</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1818.49, 1978.5)\" x=\"1818.49\" y=\"1978.5\">40</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2220.71, 1978.5)\" x=\"2220.71\" y=\"1978.5\">50</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 1872.86)\" x=\"165.504\" y=\"1872.86\">-7.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 1750.34)\" x=\"165.504\" y=\"1750.34\">-5.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 1627.82)\" x=\"165.504\" y=\"1627.82\">-2.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 1505.3)\" x=\"165.504\" y=\"1505.3\">0.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 165.504, 1382.78)\" x=\"165.504\" y=\"1382.78\">2.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(0, 1255.38, 2390.4)\" x=\"1255.38\" y=\"2390.4\">Time</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip7504)\" style=\"stroke:#0000ff; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  249.837,1487.8 290.059,1544.81 330.28,1598.37 370.502,1638.24 410.724,1685.77 450.946,1729.67 491.168,1764.38 531.389,1798.56 571.611,1832.46 611.833,1861.12 \n",
       "  652.055,1879.1 692.277,1895.25 732.499,1898.56 772.72,1897.5 812.942,1905.34 853.164,1890.46 893.386,1876.38 933.608,1880.76 973.829,1887.06 1014.05,1871.91 \n",
       "  1054.27,1848.21 1094.49,1824.11 1134.72,1792.13 1174.94,1781.08 1215.16,1743.96 1255.38,1729.96 1295.6,1691.41 1335.83,1652.57 1376.05,1629.12 1416.27,1599.45 \n",
       "  1456.49,1561.45 1496.71,1533.3 1536.93,1502.49 1577.16,1465.9 1617.38,1431.97 1657.6,1404.27 1697.82,1375.14 1738.04,1363.33 1778.27,1341.53 1818.49,1314.73 \n",
       "  1858.71,1298.72 1898.93,1280.08 1939.15,1273.34 1979.37,1266.41 2019.6,1278.59 2059.82,1301.48 2100.04,1327.07 2140.26,1353.93 2180.48,1390.55 2220.71,1425.47 \n",
       "  2260.93,1489.1 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7504)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  249.837,1487.8 290.059,1487.8 330.28,1487.8 370.502,1487.8 410.724,1487.8 450.946,1487.8 491.168,1487.8 531.389,1487.8 571.611,1487.8 611.833,1487.8 \n",
       "  652.055,1487.8 692.277,1487.8 732.499,1487.8 772.72,1487.8 812.942,1487.8 853.164,1487.8 893.386,1487.8 933.608,1487.8 973.829,1487.8 1014.05,1487.8 \n",
       "  1054.27,1487.8 1094.49,1487.8 1134.72,1487.8 1174.94,1487.8 1215.16,1487.8 1255.38,1487.8 1295.6,1487.8 1335.83,1487.8 1376.05,1487.8 1416.27,1487.8 \n",
       "  1456.49,1487.8 1496.71,1487.8 1536.93,1487.8 1577.16,1487.8 1617.38,1487.8 1657.6,1487.8 1697.82,1487.8 1738.04,1487.8 1778.27,1487.8 1818.49,1487.8 \n",
       "  1858.71,1487.8 1898.93,1487.8 1939.15,1487.8 1979.37,1487.8 2019.6,1487.8 2059.82,1487.8 2100.04,1487.8 2140.26,1487.8 2180.48,1487.8 2220.71,1487.8 \n",
       "  2260.93,1487.8 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip7501)\" points=\"\n",
       "1870.15,1451.72 2249.26,1451.72 2249.26,1330.76 1870.15,1330.76 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1870.15,1451.72 2249.26,1451.72 2249.26,1330.76 1870.15,1330.76 1870.15,1451.72 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7501)\" style=\"stroke:#0000ff; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1894.15,1391.24 2038.15,1391.24 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7501)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 2062.15, 1408.74)\" x=\"2062.15\" y=\"1408.74\">assets</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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.0\n",
    "R = zeros(4, 4)\n",
    "Rf = zeros(4, 4); Rf[1, 1] = 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.0, 0.0, 0.0, 0.0]\n",
    "C = [σ, 0.0, 0.0, 0.0]\n",
    "\n",
    "# compute solutions and simulate\n",
    "lq = QuantEcon.LQ(Q, R, A, B, C; bet = β, capT = T, rf = Rf)\n",
    "x0 = [0.0, 1, 0, 0]\n",
    "xp, up, wp = compute_sequence(lq, x0)\n",
    "\n",
    "# convert results back to assets, consumption and income\n",
    "ap = vec(xp[1, 1:end]) # assets\n",
    "c = vec(up .+ c_bar) # consumption\n",
    "time = 1:T\n",
    "income = σ * vec(wp[1, 2:end]) + m1 * time + m2 * time.^2 # income\n",
    "\n",
    "# plot results\n",
    "p1 = plot(Vector[income, ap, c, zeros(T + 1)],\n",
    "          lab = [\"non-financial income\" \"assets\" \"consumption\" \"\"],\n",
    "          color = [:orange :blue :green :black],\n",
    "          xaxis = \"Time\", layout = (2,1),\n",
    "          bottom_margin = 20mm, size = (600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\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. The\n",
    "model is solved by combining two LQ programming problems as described in\n",
    "the lecture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"600\" viewBox=\"0 0 2400 2400\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip7700\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2000\" height=\"2000\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7701\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"2400\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip7701)\" points=\"\n",
       "0,2400 2400,2400 2400,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7702\">\n",
       "    <rect x=\"480\" y=\"240\" width=\"1681\" height=\"1681\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip7701)\" points=\"\n",
       "176.123,724.503 2321.26,724.503 2321.26,47.2441 176.123,47.2441 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7703\">\n",
       "    <rect x=\"176\" y=\"47\" width=\"2146\" height=\"678\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip7703)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  202.534,724.503 202.534,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7703)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  545.537,724.503 545.537,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7703)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  888.539,724.503 888.539,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7703)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1231.54,724.503 1231.54,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7703)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1574.54,724.503 1574.54,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7703)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1917.55,724.503 1917.55,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7703)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  2260.55,724.503 2260.55,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7703)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  176.123,694.634 2321.26,694.634 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7703)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  176.123,537.28 2321.26,537.28 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7703)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  176.123,379.926 2321.26,379.926 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7703)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  176.123,222.571 2321.26,222.571 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7703)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  176.123,65.2173 2321.26,65.2173 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  176.123,724.503 2321.26,724.503 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  176.123,724.503 176.123,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  202.534,724.503 202.534,714.345 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  545.537,724.503 545.537,714.345 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  888.539,724.503 888.539,714.345 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1231.54,724.503 1231.54,714.345 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1574.54,724.503 1574.54,714.345 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1917.55,724.503 1917.55,714.345 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2260.55,724.503 2260.55,714.345 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  176.123,694.634 208.3,694.634 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  176.123,537.28 208.3,537.28 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  176.123,379.926 208.3,379.926 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  176.123,222.571 208.3,222.571 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  176.123,65.2173 208.3,65.2173 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 202.534, 778.503)\" x=\"202.534\" y=\"778.503\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 545.537, 778.503)\" x=\"545.537\" y=\"778.503\">10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 888.539, 778.503)\" x=\"888.539\" y=\"778.503\">20</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1231.54, 778.503)\" x=\"1231.54\" y=\"778.503\">30</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1574.54, 778.503)\" x=\"1574.54\" y=\"778.503\">40</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1917.55, 778.503)\" x=\"1917.55\" y=\"778.503\">50</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2260.55, 778.503)\" x=\"2260.55\" y=\"778.503\">60</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 152.123, 712.134)\" x=\"152.123\" y=\"712.134\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 152.123, 554.78)\" x=\"152.123\" y=\"554.78\">1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 152.123, 397.426)\" x=\"152.123\" y=\"397.426\">2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 152.123, 240.071)\" x=\"152.123\" y=\"240.071\">3</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 152.123, 82.7173)\" x=\"152.123\" y=\"82.7173\">4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(0, 1248.69, 1190.4)\" x=\"1248.69\" y=\"1190.4\">Time</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip7703)\" style=\"stroke:#ffa500; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  236.834,705.336 271.135,634.546 305.435,632.598 339.735,591.121 374.035,545.026 408.336,510.185 442.636,493.838 476.936,480.507 511.236,443.902 545.537,457.213 \n",
       "  579.837,417.729 614.137,398.828 648.437,360.509 682.737,312.311 717.038,313.874 751.338,292.23 785.638,275.792 819.938,278.173 854.239,234.387 888.539,240.691 \n",
       "  922.839,197.377 957.139,187.756 991.44,177.586 1025.74,168.99 1060.04,155.098 1094.34,137.591 1128.64,118.841 1162.94,132.057 1197.24,107.49 1231.54,75.867 \n",
       "  1265.84,66.4118 1300.14,115.585 1334.44,69.7055 1368.74,101.168 1403.04,101.851 1437.34,72.9153 1471.64,72.1681 1505.94,68.7258 1540.24,76.1692 1574.54,98.9373 \n",
       "  1608.84,537.28 1643.14,537.28 1677.44,537.28 1711.74,537.28 1746.04,537.28 1780.35,537.28 1814.65,537.28 1848.95,537.28 1883.25,537.28 1917.55,537.28 \n",
       "  1951.85,537.28 1986.15,537.28 2020.45,537.28 2054.75,537.28 2089.05,537.28 2123.35,537.28 2157.65,537.28 2191.95,537.28 2226.25,537.28 2260.55,537.28 \n",
       "  \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7703)\" style=\"stroke:#008000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  236.834,401.773 271.135,408.097 305.435,408.291 339.735,412.683 374.035,415.139 408.336,414.819 442.636,413.313 476.936,413.347 511.236,415.281 545.537,415.381 \n",
       "  579.837,421.32 614.137,424.728 648.437,428.711 682.737,430.064 717.038,427.054 751.338,427.503 785.638,427.571 819.938,427.974 854.239,431.673 888.539,430.96 \n",
       "  922.839,433.977 957.139,432.323 991.44,431.49 1025.74,431.264 1060.04,431.794 1094.34,432.032 1128.64,431.198 1162.94,428.902 1197.24,430.745 1231.54,429.769 \n",
       "  1265.84,424.437 1300.14,418.649 1334.44,423.48 1368.74,420.595 1403.04,424.926 1437.34,430.358 1471.64,430.649 1505.94,431.371 1540.24,431.791 1574.54,434.144 \n",
       "  1608.84,441.875 1643.14,441.875 1677.44,441.875 1711.74,441.875 1746.04,441.875 1780.35,441.875 1814.65,441.875 1848.95,441.875 1883.25,441.875 1917.55,441.875 \n",
       "  1951.85,441.875 1986.15,441.875 2020.45,441.875 2054.75,441.875 2089.05,441.875 2123.35,441.875 2157.65,441.875 2191.95,441.875 2226.25,441.875 2260.55,441.875 \n",
       "  \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip7701)\" points=\"\n",
       "1555.79,312.204 2249.26,312.204 2249.26,130.764 1555.79,130.764 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1555.79,312.204 2249.26,312.204 2249.26,130.764 1555.79,130.764 1555.79,312.204 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#ffa500; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1579.79,191.244 1723.79,191.244 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1747.79, 208.744)\" x=\"1747.79\" y=\"208.744\">non-financial income</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#008000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1579.79,251.724 1723.79,251.724 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 1747.79, 269.224)\" x=\"1747.79\" y=\"269.224\">consumption</text>\n",
       "</g>\n",
       "<polygon clip-path=\"url(#clip7701)\" points=\"\n",
       "176.123,1924.5 2321.26,1924.5 2321.26,1247.24 176.123,1247.24 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7704\">\n",
       "    <rect x=\"176\" y=\"1247\" width=\"2146\" height=\"678\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip7704)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  203.106,1924.5 203.106,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7704)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  540.391,1924.5 540.391,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7704)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  877.677,1924.5 877.677,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7704)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1214.96,1924.5 1214.96,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7704)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1552.25,1924.5 1552.25,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7704)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1889.53,1924.5 1889.53,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7704)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  2226.82,1924.5 2226.82,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7704)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  176.123,1891.47 2321.26,1891.47 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7704)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  176.123,1778.05 2321.26,1778.05 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7704)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  176.123,1664.63 2321.26,1664.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7704)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  176.123,1551.22 2321.26,1551.22 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7704)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  176.123,1437.8 2321.26,1437.8 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7704)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  176.123,1324.39 2321.26,1324.39 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  176.123,1924.5 2321.26,1924.5 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  176.123,1924.5 176.123,1247.24 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  203.106,1924.5 203.106,1914.34 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  540.391,1924.5 540.391,1914.34 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  877.677,1924.5 877.677,1914.34 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1214.96,1924.5 1214.96,1914.34 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1552.25,1924.5 1552.25,1914.34 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1889.53,1924.5 1889.53,1914.34 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2226.82,1924.5 2226.82,1914.34 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  176.123,1891.47 208.3,1891.47 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  176.123,1778.05 208.3,1778.05 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  176.123,1664.63 208.3,1664.63 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  176.123,1551.22 208.3,1551.22 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  176.123,1437.8 208.3,1437.8 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  176.123,1324.39 208.3,1324.39 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 203.106, 1978.5)\" x=\"203.106\" y=\"1978.5\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 540.391, 1978.5)\" x=\"540.391\" y=\"1978.5\">10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 877.677, 1978.5)\" x=\"877.677\" y=\"1978.5\">20</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1214.96, 1978.5)\" x=\"1214.96\" y=\"1978.5\">30</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1552.25, 1978.5)\" x=\"1552.25\" y=\"1978.5\">40</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1889.53, 1978.5)\" x=\"1889.53\" y=\"1978.5\">50</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2226.82, 1978.5)\" x=\"2226.82\" y=\"1978.5\">60</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 152.123, 1908.97)\" x=\"152.123\" y=\"1908.97\">-20</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 152.123, 1795.55)\" x=\"152.123\" y=\"1795.55\">-15</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 152.123, 1682.13)\" x=\"152.123\" y=\"1682.13\">-10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 152.123, 1568.72)\" x=\"152.123\" y=\"1568.72\">-5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 152.123, 1455.3)\" x=\"152.123\" y=\"1455.3\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 152.123, 1341.89)\" x=\"152.123\" y=\"1341.89\">5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(0, 1248.69, 2390.4)\" x=\"1248.69\" y=\"2390.4\">Time</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip7704)\" style=\"stroke:#0000ff; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  236.834,1437.8 270.563,1497.23 304.292,1537.55 338.02,1586.85 371.749,1628.46 405.477,1660.19 439.206,1686.35 472.934,1714.24 506.663,1744.77 540.391,1767.99 \n",
       "  574.12,1804.15 607.849,1831.09 641.577,1857 675.306,1873.53 709.034,1876.3 742.763,1885.55 776.491,1891.33 810.22,1895.46 843.949,1905.34 877.677,1901.56 \n",
       "  911.406,1904.49 945.134,1893.29 978.863,1881.59 1012.59,1868.81 1046.32,1855.25 1080.05,1838.36 1113.78,1816.32 1147.51,1788.31 1181.23,1767.19 1214.96,1736.94 \n",
       "  1248.69,1694.39 1282.42,1648.47 1316.15,1623.03 1349.88,1578.18 1383.61,1545.72 1417.33,1512.35 1451.06,1465.43 1484.79,1416.45 1518.52,1363.9 1552.25,1311.94 \n",
       "  1585.98,1266.41 1619.71,1271.6 1653.43,1277.04 1687.16,1282.75 1720.89,1288.75 1754.62,1295.05 1788.35,1301.67 1822.08,1308.62 1855.81,1315.91 1889.53,1323.57 \n",
       "  1923.26,1331.61 1956.99,1340.05 1990.72,1348.92 2024.45,1358.23 2058.18,1368 2091.91,1378.26 2125.63,1389.04 2159.36,1400.35 2193.09,1412.23 2226.82,1424.71 \n",
       "  2260.55,1437.81 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7704)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  236.834,1437.8 270.563,1437.8 304.292,1437.8 338.02,1437.8 371.749,1437.8 405.477,1437.8 439.206,1437.8 472.934,1437.8 506.663,1437.8 540.391,1437.8 \n",
       "  574.12,1437.8 607.849,1437.8 641.577,1437.8 675.306,1437.8 709.034,1437.8 742.763,1437.8 776.491,1437.8 810.22,1437.8 843.949,1437.8 877.677,1437.8 \n",
       "  911.406,1437.8 945.134,1437.8 978.863,1437.8 1012.59,1437.8 1046.32,1437.8 1080.05,1437.8 1113.78,1437.8 1147.51,1437.8 1181.23,1437.8 1214.96,1437.8 \n",
       "  1248.69,1437.8 1282.42,1437.8 1316.15,1437.8 1349.88,1437.8 1383.61,1437.8 1417.33,1437.8 1451.06,1437.8 1484.79,1437.8 1518.52,1437.8 1552.25,1437.8 \n",
       "  1585.98,1437.8 1619.71,1437.8 1653.43,1437.8 1687.16,1437.8 1720.89,1437.8 1754.62,1437.8 1788.35,1437.8 1822.08,1437.8 1855.81,1437.8 1889.53,1437.8 \n",
       "  1923.26,1437.8 1956.99,1437.8 1990.72,1437.8 2024.45,1437.8 2058.18,1437.8 2091.91,1437.8 2125.63,1437.8 2159.36,1437.8 2193.09,1437.8 2226.82,1437.8 \n",
       "  2260.55,1437.8 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip7701)\" points=\"\n",
       "1870.15,1451.72 2249.26,1451.72 2249.26,1330.76 1870.15,1330.76 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1870.15,1451.72 2249.26,1451.72 2249.26,1330.76 1870.15,1330.76 1870.15,1451.72 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7701)\" style=\"stroke:#0000ff; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1894.15,1391.24 2038.15,1391.24 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7701)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 2062.15, 1408.74)\" x=\"2062.15\" y=\"1408.74\">assets</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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.0\n",
    "R = zeros(4, 4)\n",
    "Rf = zeros(4, 4); Rf[1, 1] = 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.0, 0, 0, 0]\n",
    "C = zeros(4)\n",
    "\n",
    "# initialize LQ instance for retired agent\n",
    "lq_retired = QuantEcon.LQ(Q, R, A, B, C; bet = β, capT = T - K, rf = Rf)\n",
    "\n",
    "# since update_values!() changes its argument in place, we need another identical instance\n",
    "# just to get the correct value function\n",
    "lq_retired_proxy = QuantEcon.LQ(Q, R, A, B, C; bet = β, capT = T - K, rf = Rf)\n",
    "\n",
    "# iterate back to start of retirement, record final value function\n",
    "for i in 1:(T - K)\n",
    "    update_values!(lq_retired_proxy)\n",
    "end\n",
    "Rf2 = lq_retired_proxy.P\n",
    "\n",
    "# formulate LQ problem 2 (working life)\n",
    "Q = 1.0\n",
    "R = 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.0, 0, 0, 0]\n",
    "C = [σ, 0, 0, 0]\n",
    "\n",
    "# set up working life LQ instance with terminal Rf from lq_retired\n",
    "lq_working = QuantEcon.LQ(Q, R, A, B, C; bet = β, capT = K, rf = Rf2)\n",
    "\n",
    "# simulate working state / control paths\n",
    "x0 = [0.0, 1, 0, 0]\n",
    "xp_w, up_w, wp_w = compute_sequence(lq_working, x0)\n",
    "\n",
    "# simulate retirement paths (note the initial condition)\n",
    "xp_r, up_r, wp_r = compute_sequence(lq_retired, xp_w[:, end])\n",
    "\n",
    "# convert results back to assets, consumption and income\n",
    "xp = [xp_w xp_r[:, 2:end]]\n",
    "assets = vec(xp[1, :]) # assets\n",
    "\n",
    "up = [up_w up_r]\n",
    "c = vec(up .+ c_bar) # consumption\n",
    "\n",
    "time = 1:K\n",
    "income_w = σ * vec(wp_w[1, 2:K+1]) + m1 .* time + m2 .* time.^2 # income\n",
    "income_r = ones(T - K) * s\n",
    "income = [income_w; income_r]\n",
    "\n",
    "# plot results\n",
    "p2 = plot([income, assets, c, zeros(T + 1)],\n",
    "          lab = [\"non-financial income\" \"assets\" \"consumption\" \"\"],\n",
    "          color = [:orange :blue :green :black],\n",
    "          xaxis = \"Time\", layout = (2, 1),\n",
    "          bottom_margin = 20mm, size = (600, 600))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3\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": 7,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"700\" height=\"500\" viewBox=\"0 0 2800 2000\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip7900\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2000\" height=\"2000\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7901\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2800\" height=\"2000\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip7901)\" points=\"\n",
       "0,2000 2800,2000 2800,0 0,0 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7902\">\n",
       "    <rect x=\"560\" y=\"0\" width=\"1961\" height=\"1961\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polygon clip-path=\"url(#clip7901)\" points=\"\n",
       "161.394,1524.5 2721.26,1524.5 2721.26,267.717 161.394,267.717 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip7903\">\n",
       "    <rect x=\"161\" y=\"267\" width=\"2561\" height=\"1258\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip7903)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  217.743,1524.5 217.743,267.717 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7903)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1022.73,1524.5 1022.73,267.717 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7903)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1827.72,1524.5 1827.72,267.717 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7903)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  2632.71,1524.5 2632.71,267.717 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7903)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  161.394,1488.93 2721.26,1488.93 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7903)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  161.394,1155.15 2721.26,1155.15 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7903)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  161.394,821.369 2721.26,821.369 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7903)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  161.394,487.587 2721.26,487.587 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7901)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,1524.5 2721.26,1524.5 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7901)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,1524.5 161.394,267.717 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7901)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  217.743,1524.5 217.743,1505.65 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7901)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1022.73,1524.5 1022.73,1505.65 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7901)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1827.72,1524.5 1827.72,1505.65 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7901)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2632.71,1524.5 2632.71,1505.65 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7901)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,1488.93 199.792,1488.93 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7901)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,1155.15 199.792,1155.15 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7901)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,821.369 199.792,821.369 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7901)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  161.394,487.587 199.792,487.587 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7901)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 217.743, 1580.5)\" x=\"217.743\" y=\"1580.5\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7901)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1022.73, 1580.5)\" x=\"1022.73\" y=\"1580.5\">50</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7901)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1827.72, 1580.5)\" x=\"1827.72\" y=\"1580.5\">100</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7901)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2632.71, 1580.5)\" x=\"2632.71\" y=\"1580.5\">150</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7901)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 133.394, 1506.43)\" x=\"133.394\" y=\"1506.43\">2.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7901)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 133.394, 1172.65)\" x=\"133.394\" y=\"1172.65\">2.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7901)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 133.394, 838.869)\" x=\"133.394\" y=\"838.869\">3.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7901)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 133.394, 505.087)\" x=\"133.394\" y=\"505.087\">3.5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7901)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:84px; text-anchor:middle;\" transform=\"rotate(0, 1441.33, 73.2)\" x=\"1441.33\" y=\"73.2\">Dynamics with ? = 1.0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip7901)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(0, 1441.33, 1990.4)\" x=\"1441.33\" y=\"1990.4\">Time</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip7903)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  233.843,821.369 249.943,801.306 266.043,766.938 282.142,822.84 298.242,882.402 314.342,840.868 330.442,983.58 346.542,1068.18 362.641,921.374 378.741,869.229 \n",
       "  394.841,973.413 410.941,850.265 427.04,880.999 443.14,1097.03 459.24,1165.57 475.34,1420.65 491.44,1372.5 507.539,1348.77 523.639,1366.91 539.739,1347.87 \n",
       "  555.839,1343.34 571.939,1242.68 588.038,1211.75 604.138,1089.09 620.238,1033.19 636.338,995.427 652.437,996.921 668.537,908.472 684.637,843.286 700.737,866.68 \n",
       "  716.837,840.543 732.936,792.993 749.036,795.138 765.136,871.318 781.236,863.466 797.336,686.716 813.435,563.418 829.535,571.603 845.635,646.526 861.735,807.102 \n",
       "  877.834,826.419 893.934,932.597 910.034,949.16 926.134,960.407 942.234,1022.9 958.333,1007.35 974.433,848.945 990.533,924.662 1006.63,813.68 1022.73,825.887 \n",
       "  1038.83,653.275 1054.93,763.809 1071.03,904.344 1087.13,857.018 1103.23,756.198 1119.33,849.571 1135.43,751.6 1151.53,631.236 1167.63,763.897 1183.73,523.72 \n",
       "  1199.83,509.058 1215.93,527.123 1232.03,520.668 1248.13,303.286 1264.23,383.333 1280.33,605.689 1296.43,578.158 1312.53,685.84 1328.63,590.519 1344.73,563.399 \n",
       "  1360.83,577.375 1376.93,698.428 1393.03,884.483 1409.13,802.926 1425.23,822.944 1441.33,817.537 1457.43,688.325 1473.53,623.748 1489.63,517.787 1505.73,524.572 \n",
       "  1521.83,609.449 1537.93,807.667 1554.03,711.097 1570.13,829.336 1586.23,885.157 1602.32,815.502 1618.42,1000.64 1634.52,935.07 1650.62,888.621 1666.72,889.728 \n",
       "  1682.82,934.87 1698.92,762.365 1715.02,865.196 1731.12,912.996 1747.22,1103.93 1763.32,1059.24 1779.42,916.989 1795.52,872.35 1811.62,771.67 1827.72,720.952 \n",
       "  1843.82,752.93 1859.92,818.654 1876.02,985.9 1892.12,1010.75 1908.22,849.044 1924.32,691.507 1940.42,659.534 1956.52,654.335 1972.62,541.403 1988.72,419.155 \n",
       "  2004.82,396.927 2020.92,409.488 2037.02,451.375 2053.12,677.331 2069.22,695.674 2085.32,851.93 2101.42,779.258 2117.52,749.67 2133.62,854.541 2149.72,893.71 \n",
       "  2165.82,1010.63 2181.92,975.509 2198.02,880.691 2214.12,952.045 2230.22,1009.6 2246.32,1184.33 2262.42,1161.65 2278.52,1299.3 2294.62,1212.54 2310.72,1220.49 \n",
       "  2326.82,1226.81 2342.91,1012.98 2359.01,1001.43 2375.11,826.644 2391.21,949.837 2407.31,878.45 2423.41,850.724 2439.51,920.962 2455.61,935.966 2471.71,1088.59 \n",
       "  2487.81,1035.5 2503.91,947.528 2520.01,914.588 2536.11,698.889 2552.21,689.479 2568.31,459.551 2584.41,499.672 2600.51,589.615 2616.61,601.665 2632.71,596.69 \n",
       "  2648.81,654.995 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7903)\" style=\"stroke:#0000ff; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  233.843,1488.93 249.943,1166.59 266.043,991.946 282.142,888.01 298.242,856.415 314.342,863.68 330.442,850.977 346.542,900.965 362.641,960.343 378.741,932.87 \n",
       "  394.841,897.998 410.941,921.252 427.04,884.474 443.14,877.635 459.24,959.71 475.34,1029.32 491.44,1166.4 507.539,1218.21 523.639,1235.6 539.739,1251.78 \n",
       "  555.839,1252.61 571.939,1251.24 588.038,1210.64 604.138,1177.38 620.238,1111.58 636.338,1055.39 652.437,1011.37 668.537,989.198 684.637,942.679 700.737,892.789 \n",
       "  716.837,876.26 732.936,857.354 749.036,828.733 765.136,814.781 781.236,837.757 797.336,846.527 813.435,781.016 829.535,698.274 845.635,658.729 861.735,667.971 \n",
       "  877.834,736.387 893.934,779.423 910.034,843.757 926.134,883.591 942.234,908.648 958.333,946.374 974.433,959.72 990.533,903.844 1006.63,904.955 1022.73,861.547 \n",
       "  1038.83,843.937 1054.93,766.424 1071.03,770.144 1087.13,827.762 1103.23,838.803 1119.33,804.557 1135.43,823.852 1151.53,795.003 1167.63,732.384 1183.73,752.575 \n",
       "  1199.83,667.834 1215.93,618.2 1232.03,599.692 1248.13,587.563 1264.23,495.141 1280.33,479.068 1296.43,558.877 1312.53,589.239 1328.63,647.615 1344.73,640.028 \n",
       "  1360.83,625.356 1376.93,623.308 1393.03,670.222 1409.13,768.217 1425.23,786.573 1441.33,803.999 1457.43,810.868 1473.53,763.213 1489.63,712.976 1505.73,645.005 \n",
       "  1521.83,612.543 1537.93,629.392 1554.03,716.66 1570.13,723.519 1586.23,773.924 1602.32,822.113 1618.42,819.428 1634.52,891.412 1650.62,902.651 1666.72,890.055 \n",
       "  1682.82,883.98 1698.92,898.729 1715.02,837.991 1731.12,847.334 1747.22,871.108 1763.32,959.072 1779.42,986.849 1795.52,944.84 1811.62,905.424 1827.72,845.141 \n",
       "  1843.82,793.867 1859.92,780.024 1876.02,798.912 1892.12,874.96 1908.22,924.135 1924.32,885.481 1940.42,803.059 1956.52,747.765 1972.62,717.109 1988.72,656.501 \n",
       "  2004.82,576.71 2020.92,526.639 2037.02,505.723 2053.12,511.506 2069.22,604.044 2085.32,659.169 2101.42,749.6 2117.52,767.566 2133.62,765.131 2149.72,805.432 \n",
       "  2165.82,841.796 2181.92,906.937 2198.02,926.706 2214.12,899.352 2230.22,913.484 2246.32,943.603 2262.42,1028.42 2278.52,1063.3 2294.62,1135.88 2310.72,1139.04 \n",
       "  2326.82,1143.82 2342.91,1148.8 2359.01,1066.63 2375.11,1019.56 2391.21,925.953 2407.31,926.365 2423.41,898.288 2439.51,872.78 2455.61,887.425 2471.71,900.944 \n",
       "  2487.81,968.421 2503.91,982.277 2520.01,954.578 2536.11,927.2 2552.21,827.559 2568.31,772.301 2584.41,652.605 2600.51,606.605 2616.61,618.461 2632.71,629.368 \n",
       "  2648.81,633.037 \n",
       "  \"/>\n",
       "<polygon clip-path=\"url(#clip7901)\" points=\"\n",
       "2247.56,556.677 2637.26,556.677 2637.26,375.237 2247.56,375.237 \n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip7901)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2247.56,556.677 2637.26,556.677 2637.26,375.237 2247.56,375.237 2247.56,556.677 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip7901)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2275.56,435.717 2443.56,435.717 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7901)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 2471.56, 453.217)\" x=\"2471.56\" y=\"453.217\">q bar</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip7901)\" style=\"stroke:#0000ff; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2275.56,496.197 2443.56,496.197 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip7901)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 2471.56, 513.697)\" x=\"2471.56\" y=\"513.697\">q</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# model parameters\n",
    "a0 = 5.0\n",
    "a1 = 0.5\n",
    "σ = 0.15\n",
    "ρ = 0.9\n",
    "γ = 1.0\n",
    "β = 0.95\n",
    "c = 2.0\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; -a1 a1 0; 0 0 0]\n",
    "A = [ρ 0 m0 * (1 - ρ); 0 1 0; 0 0 1]\n",
    "\n",
    "B = [0.0, 1, 0]\n",
    "C = [m1 * σ, 0, 0]\n",
    "\n",
    "lq = QuantEcon.LQ(Q, R, A, B, C; bet = β)\n",
    "\n",
    "# simulate state / control paths\n",
    "x0 = [m0, 2, 1]\n",
    "xp, up, wp = compute_sequence(lq, x0, 150)\n",
    "q_bar = vec(xp[1, :])\n",
    "q = vec(xp[2, :])\n",
    "\n",
    "# plot simulation results\n",
    "p3 = plot(eachindex(q), [q_bar q],\n",
    "          lab = [\"q bar\" \"q\"],\n",
    "          color = [:black :blue],\n",
    "          xaxis = \"Time\", title = \"Dynamics with γ = $γ\",\n",
    "          bottom_margin = 20mm, top_margin = 10mm,\n",
    "          size = (700, 500))"
   ]
  }
 ],
 "metadata": {
  "filename": "lqcontrol.rst",
  "kernelspec": {
   "display_name": "Julia",
   "language": "julia",
   "name": "julia-1.0"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.0.2"
  },
  "title": "LQ Dynamic Programming Problems"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}