lu_tricks.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='lu-tricks'></a>\n",
    "<div id=\"qe-notebook-header\" align=\"right\" style=\"text-align:right;\">\n",
    "        <a href=\"https://quantecon.org/\" title=\"quantecon.org\">\n",
    "                <img style=\"width:250px;display:inline;\" width=\"250px\" src=\"https://assets.quantecon.org/img/qe-menubar-logo.svg\" alt=\"QuantEcon\">\n",
    "        </a>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Classical Control with Linear Algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Classical Control with Linear Algebra](#Classical-Control-with-Linear-Algebra)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [A Control Problem](#A-Control-Problem)  \n",
    "  - [Finite Horizon Theory](#Finite-Horizon-Theory)  \n",
    "  - [The Infinite Horizon Limit](#The-Infinite-Horizon-Limit)  \n",
    "  - [Undiscounted Problems](#Undiscounted-Problems)  \n",
    "  - [Implementation](#Implementation)  \n",
    "  - [Exercises](#Exercises)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "In an earlier lecture [Linear Quadratic Dynamic Programming Problems](https://python.quantecon.org/lqcontrol.html), we have studied how to solve a special\n",
    "class of dynamic optimization and prediction problems by applying the method of dynamic programming. In this class of problems\n",
    "\n",
    "> - the objective function is **quadratic** in **states** and **controls**.  \n",
    "- the one-step transition function is **linear**.  \n",
    "- shocks are IID Gaussian or martingale differences.  \n",
    "\n",
    "\n",
    "\n",
    "In this lecture and a companion lecture [Classical Filtering with Linear Algebra](https://python.quantecon.org/classical_filtering.html), we study the classical theory of linear-quadratic (LQ) optimal control problems.\n",
    "\n",
    "The classical approach does not use the two closely related methods – dynamic programming and  Kalman filtering – that we describe in other lectures, namely, [Linear Quadratic Dynamic Programming Problems](https://python.quantecon.org/lqcontrol.html) and [A First Look at the Kalman Filter](https://python.quantecon.org/kalman.html).\n",
    "\n",
    "Instead, they use either\n",
    "\n",
    "> - $ z $-transform and lag operator methods, or  \n",
    "- matrix decompositions applied to linear systems of first-order conditions for optimum problems.  \n",
    "\n",
    "\n",
    "\n",
    "In this lecture and the sequel [Classical Filtering with Linear Algebra](https://python.quantecon.org/classical_filtering.html), we mostly rely on  elementary linear algebra.\n",
    "\n",
    "The main tool from linear algebra we’ll put to work here is [LU decomposition](https://en.wikipedia.org/wiki/LU_decomposition).\n",
    "\n",
    "We’ll begin with discrete horizon problems.\n",
    "\n",
    "Then we’ll view infinite horizon problems as appropriate limits of these finite horizon problems.\n",
    "\n",
    "Later, we will examine the close connection between LQ control and least-squares prediction and filtering problems.\n",
    "\n",
    "These classes of problems are connected in the sense that to solve each, essentially the same mathematics is used.\n",
    "\n",
    "Let’s start with some standard imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "Useful references include [[Whi63]](https://python.quantecon.org/zreferences.html#whittle1963), [[HS80]](https://python.quantecon.org/zreferences.html#hansar1980), [[Orf88]](https://python.quantecon.org/zreferences.html#orfanidisoptimum1988), [[AP91]](https://python.quantecon.org/zreferences.html#athanasios1991), and [[Mut60]](https://python.quantecon.org/zreferences.html#muth1960)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Control Problem\n",
    "\n",
    "Let $ L $ be the **lag operator**, so that, for sequence $ \\{x_t\\} $ we have $ L x_t = x_{t-1} $.\n",
    "\n",
    "More generally, let $ L^k x_t = x_{t-k} $ with $ L^0 x_t = x_t $ and\n",
    "\n",
    "$$\n",
    "d(L) = d_0 + d_1 L+ \\ldots + d_m L^m\n",
    "$$\n",
    "\n",
    "where $ d_0, d_1, \\ldots, d_m $ is a given scalar sequence.\n",
    "\n",
    "Consider the discrete-time control problem\n",
    "\n",
    "\n",
    "<a id='equation-oneone'></a>\n",
    "$$\n",
    "\\max_{\\{y_t\\}}\n",
    "\\lim_{N \\to \\infty} \\sum^N_{t=0} \\beta^t\\,\n",
    "\\left\\{\n",
    "     a_t y_t - {1 \\over 2}\\, hy^2_t - {1 \\over 2} \\,\n",
    "         \\left[ d(L)y_t \\right]^2\n",
    "\\right\\}, \\tag{1}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "- $ h $ is a positive parameter and $ \\beta \\in (0,1) $ is a discount factor.  \n",
    "- $ \\{a_t\\}_{t \\geq 0} $ is a sequence of exponential order less than $ \\beta^{-1/2} $, by which we mean $ \\lim_{t \\rightarrow \\infty} \\beta^{\\frac{t}{2}} a_t = 0 $.  \n",
    "\n",
    "\n",
    "Maximization in [(1)](#equation-oneone) is subject to  initial conditions for $ y_{-1}, y_{-2} \\ldots, y_{-m} $.\n",
    "\n",
    "Maximization is over infinite sequences $ \\{y_t\\}_{t \\geq 0} $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "The formulation of the LQ problem given above is broad enough to encompass\n",
    "many useful models.\n",
    "\n",
    "As a simple illustration, recall that in [LQ Control: Foundations](https://python.quantecon.org/lqcontrol.html) we consider a monopolist facing stochastic demand\n",
    "shocks and adjustment costs.\n",
    "\n",
    "Let’s consider a deterministic version of this problem, where the monopolist\n",
    "maximizes the discounted sum\n",
    "\n",
    "$$\n",
    "\\sum_{t=0}^{\\infty} \\beta^t \\pi_t\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "\\pi_t = p_t q_t - c q_t - \\gamma (q_{t+1} - q_t)^2\n",
    "\\quad \\text{with} \\quad\n",
    "p_t = \\alpha_0 - \\alpha_1 q_t + d_t\n",
    "$$\n",
    "\n",
    "In this expression, $ q_t $ is output, $ c $ is average cost of production, and $ d_t $ is a demand shock.\n",
    "\n",
    "The term $ \\gamma (q_{t+1} - q_t)^2 $ represents adjustment costs.\n",
    "\n",
    "You will be able to confirm that the objective function can be rewritten as [(1)](#equation-oneone) when\n",
    "\n",
    "- $ a_t := \\alpha_0 + d_t - c $  \n",
    "- $ h := 2 \\alpha_1 $  \n",
    "- $ d(L) := \\sqrt{2 \\gamma}(I - L) $  \n",
    "\n",
    "\n",
    "Further examples of this problem for factor demand, economic growth, and government policy problems are given in ch. IX of [[Sar87]](https://python.quantecon.org/zreferences.html#sargent1987)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finite Horizon Theory\n",
    "\n",
    "We first study a finite $ N $ version of the problem.\n",
    "\n",
    "Later we will study an infinite horizon problem solution as a limiting version of a finite horizon problem.\n",
    "\n",
    "(This will require being careful because the limits as $ N \\to \\infty $ of the necessary and sufficient conditions for maximizing finite $ N $ versions of [(1)](#equation-oneone)\n",
    "are not sufficient for maximizing [(1)](#equation-oneone))\n",
    "\n",
    "We begin by\n",
    "\n",
    "1. fixing $ N > m $,  \n",
    "1. differentiating the finite version of [(1)](#equation-oneone) with respect to $ y_0, y_1, \\ldots, y_N $, and  \n",
    "1. setting these derivatives to zero.  \n",
    "\n",
    "\n",
    "For $ t=0, \\ldots, N-m $ these first-order necessary conditions are the\n",
    "*Euler equations*.\n",
    "\n",
    "For $ t = N-m + 1, \\ldots, N $, the first-order conditions are a set of\n",
    "*terminal conditions*.\n",
    "\n",
    "Consider the term\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "J\n",
    "& = \\sum^N_{t=0} \\beta^t [d(L) y_t] [d(L) y_t]\n",
    "\\\\\n",
    "& = \\sum^N_{t=0}\n",
    "    \\beta^t \\, (d_0 \\, y_t + d_1 \\, y_{t-1} + \\cdots + d_m \\, y_{t-m}) \\,\n",
    "               (d_0 \\, y_t + d_1 \\, y_{t-1} + \\cdots  + d_m\\, y_{t-m})\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Differentiating $ J $ with respect to $ y_t $ for\n",
    "$ t=0,\\ 1,\\ \\ldots,\\ N-m $ gives\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "{\\partial {J} \\over \\partial y_t}\n",
    "   & = 2 \\beta^t \\, d_0 \\, d(L)y_t +\n",
    "       2 \\beta^{t+1} \\, d_1\\, d(L)y_{t+1} + \\cdots +\n",
    "       2 \\beta^{t+m}\\, d_m\\, d(L) y_{t+m} \\\\\n",
    "   & = 2\\beta^t\\, \\bigl(d_0 + d_1 \\, \\beta L^{-1} + d_2 \\, \\beta^2\\, L^{-2} +\n",
    "       \\cdots + d_m \\, \\beta^m \\, L^{-m}\\bigr)\\, d (L) y_t\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "We can write this more succinctly as\n",
    "\n",
    "\n",
    "<a id='equation-onetwo'></a>\n",
    "$$\n",
    "{\\partial {J} \\over \\partial y_t}\n",
    "    = 2 \\beta^t \\, d(\\beta L^{-1}) \\, d (L) y_t \\tag{2}\n",
    "$$\n",
    "\n",
    "Differentiating $ J $ with respect to $ y_t $ for $ t = N-m + 1, \\ldots, N $ gives\n",
    "\n",
    "\n",
    "<a id='equation-onethree'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    " {\\partial J \\over \\partial y_N}\n",
    " &= 2 \\beta^N\\, d_0 \\, d(L) y_N \\cr\n",
    "   {\\partial J \\over \\partial y_{N-1}}\n",
    " &= 2\\beta^{N-1} \\,\\bigl[d_0 + \\beta \\,\n",
    "   d_1\\, L^{-1}\\bigr] \\, d(L)y_{N-1} \\cr\n",
    "   \\vdots\n",
    " & \\quad \\quad \\vdots \\cr\n",
    "   {\\partial {J} \\over \\partial y_{N-m+1}}\n",
    " &= 2 \\beta^{N-m+1}\\,\\bigl[d_0 + \\beta\n",
    "   L^{-1} \\,d_1 + \\cdots + \\beta^{m-1}\\, L^{-m+1}\\, d_{m-1}\\bigr]  d(L)y_{N-m+1}\n",
    "\\end{aligned} \\tag{3}\n",
    "$$\n",
    "\n",
    "With these preliminaries under our belts, we are ready to differentiate [(1)](#equation-oneone).\n",
    "\n",
    "Differentiating [(1)](#equation-oneone) with respect to $ y_t $ for $ t=0, \\ldots, N-m $ gives the Euler equations\n",
    "\n",
    "\n",
    "<a id='equation-onefour'></a>\n",
    "$$\n",
    "\\bigl[h+d\\,(\\beta L^{-1})\\,d(L)\\bigr] y_t = a_t,\n",
    "\\quad t=0,\\, 1,\\, \\ldots, N-m \\tag{4}\n",
    "$$\n",
    "\n",
    "The system of equations [(4)](#equation-onefour) forms  a  $ 2 \\times m $ order linear *difference\n",
    "equation* that must hold for the values of $ t $ indicated.\n",
    "\n",
    "Differentiating [(1)](#equation-oneone) with respect to $ y_t $ for $ t = N-m + 1, \\ldots, N $ gives the terminal conditions\n",
    "\n",
    "\n",
    "<a id='equation-onefive'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\beta^N  (a_N - hy_N - d_0\\,d(L)y_N)\n",
    "&= 0  \\cr\n",
    "  \\beta^{N-1} \\left(a_{N-1}-hy_{N-1}-\\Bigl(d_0 + \\beta \\, d_1\\,\n",
    "L^{-1}\\Bigr)\\, d(L)\\, y_{N-1}\\right)\n",
    "& = 0 \\cr\n",
    " \\vdots & \\vdots\\cr\n",
    "\\beta^{N-m+1} \\biggl(a_{N-m+1} - h y_{N-m+1} -(d_0+\\beta L^{-1}\n",
    "d_1+\\cdots\\  +\\beta^{m-1} L^{-m+1} d_{m-1}) d(L) y_{N-m+1}\\biggr)\n",
    "& = 0\n",
    "\\end{aligned} \\tag{5}\n",
    "$$\n",
    "\n",
    "In the finite $ N $ problem, we want simultaneously to solve [(4)](#equation-onefour) subject to the $ m $ initial conditions\n",
    "$ y_{-1}, \\ldots, y_{-m} $ and the $ m $ terminal conditions\n",
    "[(5)](#equation-onefive).\n",
    "\n",
    "These conditions uniquely pin down the solution of the finite $ N $ problem.\n",
    "\n",
    "That is, for the finite $ N $ problem,\n",
    "conditions [(4)](#equation-onefour) and [(5)](#equation-onefive) are necessary and sufficient for a maximum,\n",
    "by concavity of the objective function.\n",
    "\n",
    "Next, we describe how to obtain the solution using matrix methods.\n",
    "\n",
    "\n",
    "<a id='fdlq'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix Methods\n",
    "\n",
    "Let’s look at how linear algebra can be used to tackle and shed light on the finite horizon LQ control problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A Single Lag Term\n",
    "\n",
    "Let’s begin with the special case in which $ m=1 $.\n",
    "\n",
    "We want to solve the system of $ N+1 $ linear equations\n",
    "\n",
    "\n",
    "<a id='equation-oneff'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\bigl[h & + d\\, (\\beta L^{-1})\\, d\\, (L) ] y_t = a_t, \\quad\n",
    "t = 0,\\ 1,\\ \\ldots,\\, N-1\\cr\n",
    "\\beta^N & \\bigl[a_N-h\\, y_N-d_0\\, d\\, (L) y_N\\bigr] = 0\n",
    "\\end{aligned} \\tag{6}\n",
    "$$\n",
    "\n",
    "where $ d(L) = d_0 + d_1 L $.\n",
    "\n",
    "These equations are to be solved for\n",
    "$ y_0, y_1, \\ldots, y_N $ as functions of\n",
    "$ a_0, a_1, \\ldots,  a_N $ and $ y_{-1} $.\n",
    "\n",
    "Let\n",
    "\n",
    "$$\n",
    "\\phi (L)\n",
    "= \\phi_0 + \\phi_1 L + \\beta \\phi_1 L^{-1}\n",
    "= h + d (\\beta L^{-1}) d(L)\n",
    "= (h + d_0^2 + d_1^2) + d_1 d_0 L+ d_1 d_0 \\beta L^{-1}\n",
    "$$\n",
    "\n",
    "Then we can represent [(6)](#equation-oneff) as the matrix equation\n",
    "\n",
    "\n",
    "<a id='equation-onefourfive'></a>\n",
    "$$\n",
    "\\left[\n",
    "    \\begin{matrix}\n",
    "        (\\phi_0-d_1^2) & \\phi_1 & 0 & 0 & \\ldots & \\ldots & 0 \\cr\n",
    "        \\beta \\phi_1 & \\phi_0 & \\phi_1 & 0 & \\ldots & \\dots & 0 \\cr\n",
    "        0 & \\beta \\phi_1 & \\phi_0 & \\phi_1 & \\ldots & \\ldots & 0 \\cr\n",
    "        \\vdots &\\vdots & \\vdots & \\ddots & \\vdots & \\vdots & \\vdots \\cr\n",
    "        0 & \\ldots & \\ldots & \\ldots & \\beta \\phi_1 & \\phi_0 &\\phi_1 \\cr\n",
    "        0 & \\ldots & \\ldots & \\ldots & 0 & \\beta \\phi_1 & \\phi_0\n",
    "    \\end{matrix}\n",
    "\\right]\n",
    "\\left [\n",
    "    \\begin{matrix}\n",
    "        y_N \\cr y_{N-1} \\cr y_{N-2} \\cr \\vdots \\cr\n",
    "        y_1 \\cr y_0\n",
    "    \\end{matrix}\n",
    "\\right ] =\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "    a_N \\cr a_{N-1} \\cr a_{N-2} \\cr \\vdots \\cr a_1 \\cr\n",
    "    a_0 - \\phi_1 y_{-1}\n",
    "\\end{matrix}\n",
    "\\right] \\tag{7}\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "\n",
    "<a id='equation-onefoursix'></a>\n",
    "$$\n",
    "W\\bar y = \\bar a \\tag{8}\n",
    "$$\n",
    "\n",
    "Notice how we have chosen to arrange the $ y_t $’s in reverse\n",
    "time order.\n",
    "\n",
    "The matrix $ W $ on the left side of [(7)](#equation-onefourfive) is “almost” a\n",
    "[Toeplitz matrix](https://en.wikipedia.org/wiki/Toeplitz_matrix) (where each\n",
    "descending diagonal is constant).\n",
    "\n",
    "There are two sources of deviation from the  form  of a Toeplitz matrix\n",
    "\n",
    "1. The first element differs from the remaining diagonal elements, reflecting the terminal condition.  \n",
    "1. The sub-diagonal elements equal $ \\beta $ time the super-diagonal elements.  \n",
    "\n",
    "\n",
    "The solution of [(8)](#equation-onefoursix) can be expressed in the form\n",
    "\n",
    "\n",
    "<a id='equation-onefourseven'></a>\n",
    "$$\n",
    "\\bar y = W^{-1} \\bar a \\tag{9}\n",
    "$$\n",
    "\n",
    "which represents each element $ y_t $ of $ \\bar y $ as a function of the entire vector $ \\bar a $.\n",
    "\n",
    "That is, $ y_t $ is a function of past, present, and future values of $ a $’s, as well as of the initial condition $ y_{-1} $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### An Alternative Representation\n",
    "\n",
    "An alternative way to express the solution to [(7)](#equation-onefourfive) or\n",
    "[(8)](#equation-onefoursix) is in so-called **feedback-feedforward** form.\n",
    "\n",
    "The idea here is to find a solution expressing $ y_t $ as a function of *past* $ y $’s and *current* and *future* $ a $’s.\n",
    "\n",
    "To achieve this solution, one can use an [LU decomposition](https://en.wikipedia.org/wiki/LU_decomposition) of $ W $.\n",
    "\n",
    "There always exists a decomposition of $ W $ of the form $ W= LU $\n",
    "where\n",
    "\n",
    "- $ L $ is an $ (N+1) \\times (N+1) $ lower triangular matrix.  \n",
    "- $ U $ is an $ (N+1) \\times (N+1) $ upper triangular matrix.  \n",
    "\n",
    "\n",
    "The factorization can be normalized so that the diagonal elements of $ U $ are unity.\n",
    "\n",
    "Using the LU representation in [(9)](#equation-onefourseven), we obtain\n",
    "\n",
    "\n",
    "<a id='equation-onefournine'></a>\n",
    "$$\n",
    "U \\bar y = L^{-1} \\bar a \\tag{10}\n",
    "$$\n",
    "\n",
    "Since $ L^{-1} $ is lower triangular, this representation expresses\n",
    "$ y_t $ as a function of\n",
    "\n",
    "- lagged $ y $’s (via the term $ U \\bar y $), and  \n",
    "- current and future $ a $’s (via the term $ L^{-1} \\bar a $)  \n",
    "\n",
    "\n",
    "Because there are zeros everywhere in the matrix\n",
    "on the left of [(7)](#equation-onefourfive) except on the diagonal, super-diagonal, and\n",
    "sub-diagonal, the $ LU $ decomposition takes\n",
    "\n",
    "- $ L $ to be zero except in the diagonal  and the leading sub-diagonal.  \n",
    "- $ U $ to be zero except on the diagonal and the super-diagonal.  \n",
    "\n",
    "\n",
    "Thus, [(10)](#equation-onefournine) has the form\n",
    "\n",
    "$$\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "    1& U_{12} & 0 & 0 & \\ldots & 0 & 0 \\cr\n",
    "    0 & 1 & U_{23} & 0 & \\ldots & 0 & 0 \\cr\n",
    "    0 & 0 & 1 & U_{34} & \\ldots & 0 & 0 \\cr\n",
    "    0 & 0 & 0 & 1 & \\ldots & 0 & 0\\cr\n",
    "    \\vdots & \\vdots & \\vdots & \\vdots & \\ddots & \\vdots & \\vdots\\cr\n",
    "    0 & 0 & 0 & 0 & \\ldots & 1 & U_{N,N+1} \\cr\n",
    "    0 & 0 & 0 & 0 & \\ldots & 0 & 1\n",
    "\\end{matrix}\n",
    "\\right] \\ \\ \\\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "    y_N \\cr y_{N-1} \\cr y_{N-2} \\cr y_{N-3} \\cr \\vdots \\cr y_1 \\cr y_0\n",
    "\\end{matrix}\n",
    "\\right] =\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\quad\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "    L^{-1}_{11} & 0 & 0 & \\ldots & 0 \\cr\n",
    "    L^{-1}_{21} & L^{-1}_{22} & 0 & \\ldots & 0 \\cr\n",
    "    L^{-1}_{31} & L^{-1}_{32} & L^{-1}_{33}& \\ldots & 0 \\cr\n",
    "    \\vdots & \\vdots & \\vdots & \\ddots & \\vdots\\cr\n",
    "    L^{-1}_{N,1} & L^{-1}_{N,2} & L^{-1}_{N,3} & \\ldots & 0 \\cr\n",
    "    L^{-1}_{N+1,1} & L^{-1}_{N+1,2} & L^{-1}_{N+1,3} & \\ldots &\n",
    "    L^{-1}_{N+1\\, N+1}\n",
    "\\end{matrix}\n",
    "\\right]\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "    a_N \\cr a_{N-1} \\cr a_{N-2} \\cr \\vdots \\cr a_1 \\cr a_0 -\n",
    "    \\phi_1 y_{-1}\n",
    "\\end{matrix}\n",
    "\\right ]\n",
    "$$\n",
    "\n",
    "where $ L^{-1}_{ij} $ is the $ (i,j) $ element of $ L^{-1} $ and $ U_{ij} $ is the $ (i,j) $ element of $ U $.\n",
    "\n",
    "Note how the left side for a given $ t $ involves  $ y_t $ and one lagged value $ y_{t-1} $ while the right side involves all future values of the forcing process $ a_t, a_{t+1}, \\ldots, a_N $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Additional Lag Terms\n",
    "\n",
    "We briefly indicate how this approach extends to the problem with\n",
    "$ m > 1 $.\n",
    "\n",
    "Assume that $ \\beta = 1 $ and  let $ D_{m+1} $ be the\n",
    "$ (m+1) \\times (m+1) $ symmetric matrix whose elements are\n",
    "determined from the following formula:\n",
    "\n",
    "$$\n",
    "D_{jk} = d_0 d_{k-j} + d_1 d_{k-j+1} + \\ldots + d_{j-1} d_{k-1}, \\qquad k\n",
    "\\geq j\n",
    "$$\n",
    "\n",
    "Let $ I_{m+1} $ be the $ (m+1) \\times (m+1) $ identity matrix.\n",
    "\n",
    "Let $ \\phi_j $ be the coefficients in the expansion $ \\phi (L) = h + d (L^{-1}) d (L) $.\n",
    "\n",
    "Then the first order conditions [(4)](#equation-onefour) and [(5)](#equation-onefive) can be expressed as:\n",
    "\n",
    "$$\n",
    "(D_{m+1} + hI_{m+1})\\ \\\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "    y_N \\cr y_{N-1} \\cr \\vdots \\cr y_{N-m}\n",
    "\\end{matrix}\n",
    "\\right]\\\n",
    "= \\ \\left[\n",
    "\\begin{matrix}\n",
    "     a_N \\cr a_{N-1} \\cr \\vdots \\cr a_{N-m}\n",
    " \\end{matrix}\n",
    "\\right] + M\\\n",
    "\\left[\n",
    "    \\begin{matrix}\n",
    "        y_{N-m+1}\\cr y_{N-m-2}\\cr \\vdots\\cr y_{N-2m}\n",
    "    \\end{matrix}\n",
    "\\right]\n",
    "$$\n",
    "\n",
    "where $ M $ is $ (m+1)\\times m $ and\n",
    "\n",
    "$$\n",
    "M_{ij} = \\begin{cases}\n",
    "D_{i-j,\\,m+1} \\textrm{ for } i>j  \\\\\n",
    " 0  \\textrm{ for }  i\\leq j\\end{cases}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\phi_m y_{N-1} &+ \\phi_{m-1} y_{N-2} + \\ldots + \\phi_0 y_{N-m-1} +\n",
    "\\phi_1 y_{N-m-2} +\\cr\n",
    "&\\hskip.75in \\ldots + \\phi_m y_{N-2m-1} = a_{N-m-1} \\cr\n",
    "\\phi_m y_{N-2} &+ \\phi_{m-1} y_{N-3} + \\ldots + \\phi_0 y_{N-m-2} + \\phi_1\n",
    "y_{N-m-3} +\\cr\n",
    "&\\hskip.75in \\ldots + \\phi_m y_{N-2m-2} = a_{N-m-2} \\cr\n",
    "&\\qquad \\vdots \\cr\n",
    "\\phi_m y_{m+1} &+ \\phi_{m-1} y_m + + \\ldots + \\phi_0 y_1 + \\phi_1 y_0 +\n",
    "\\phi_m y_{-m+1} = a_1 \\cr\n",
    "\\phi_m y_m + \\phi_{m-1}& y_{m-1} + \\phi_{m-2} + \\ldots + \\phi_0 y_0 + \\phi_1\n",
    "y_{-1} + \\ldots + \\phi_m y_{-m} = a_0\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "As before, we can express this equation as $ W \\bar y = \\bar a $.\n",
    "\n",
    "The matrix on the left of this equation is “almost” Toeplitz, the\n",
    "exception being the leading $ m \\times m $ submatrix in the upper\n",
    "left-hand corner.\n",
    "\n",
    "We can represent the solution in feedback-feedforward form by obtaining a decomposition $ LU = W $, and obtain\n",
    "\n",
    "\n",
    "<a id='equation-onefivetwo'></a>\n",
    "$$\n",
    "U \\bar y = L^{-1} \\bar a \\tag{11}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{aligned} \\sum^t_{j=0}\\, U_{-t+N+1,\\,-t+N+j+1}\\,y_{t-j} &= \\sum^{N-t}_{j=0}\\,\n",
    "L_{-t+N+1,\\, -t+N+1-j}\\, \\bar a_{t+j}\\ ,\\cr\n",
    "&\\qquad t=0,1,\\ldots, N\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "where $ L^{-1}_{t,s} $ is the element in the $ (t,s) $ position\n",
    "of $ L $, and similarly for $ U $.\n",
    "\n",
    "The left side of equation [(11)](#equation-onefivetwo) is the “feedback” part of the optimal\n",
    "control law for $ y_t $, while the right-hand side is the “feedforward” part.\n",
    "\n",
    "We note that there is a different control law for each $ t $.\n",
    "\n",
    "Thus, in the finite horizon case, the optimal control law is time-dependent.\n",
    "\n",
    "It is natural to suspect that as $ N \\rightarrow\\infty $, [(11)](#equation-onefivetwo)\n",
    "becomes equivalent to the solution of our infinite horizon problem,\n",
    "which below we shall show can be expressed as\n",
    "\n",
    "$$\n",
    "c(L) y_t = c (\\beta L^{-1})^{-1} a_t\\ ,\n",
    "$$\n",
    "\n",
    "so that as $ N \\rightarrow \\infty $ we expect that for each fixed\n",
    "$ t, U^{-1}_{t, t-j}\n",
    "\\rightarrow c_j $ and $ L_{t,t+j} $ approaches the coefficient on\n",
    "$ L^{-j} $ in the expansion of $ c(\\beta L^{-1})^{-1} $.\n",
    "\n",
    "This suspicion is true under general conditions that we shall study later.\n",
    "\n",
    "For now, we note that by creating the matrix $ W $ for large\n",
    "$ N $ and factoring it into the $ LU $ form, good approximations\n",
    "to $ c(L) $ and $ c(\\beta L^{-1})^{-1} $ can be obtained."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Infinite Horizon Limit\n",
    "\n",
    "For the infinite horizon problem, we propose to discover first-order\n",
    "necessary conditions by taking the limits of [(4)](#equation-onefour) and [(5)](#equation-onefive) as\n",
    "$ N \\to \\infty $.\n",
    "\n",
    "This approach is valid, and the limits of [(4)](#equation-onefour) and [(5)](#equation-onefive) as $ N $ approaches infinity are first-order necessary conditions for a maximum.\n",
    "\n",
    "However, for the infinite horizon problem with $ \\beta < 1 $, the limits of [(4)](#equation-onefour) and [(5)](#equation-onefive) are, in general, not sufficient for a maximum.\n",
    "\n",
    "That is, the limits of [(5)](#equation-onefive) do not provide enough information uniquely to determine the solution of the Euler equation [(4)](#equation-onefour) that maximizes [(1)](#equation-oneone).\n",
    "\n",
    "As we shall see below, a side condition on the path of $ y_t $ that together with [(4)](#equation-onefour) is sufficient for an optimum is\n",
    "\n",
    "\n",
    "<a id='equation-onesix'></a>\n",
    "$$\n",
    "\\sum^\\infty_{t=0}\\ \\beta^t\\, hy^2_t < \\infty \\tag{12}\n",
    "$$\n",
    "\n",
    "All paths that satisfy the Euler equations, except the one that we shall\n",
    "select below, violate this condition and, therefore, evidently lead to\n",
    "(much) lower values of [(1)](#equation-oneone) than does the\n",
    "optimal path selected by the solution procedure below.\n",
    "\n",
    "Consider the *characteristic equation* associated with the Euler equation\n",
    "\n",
    "\n",
    "<a id='equation-oneseven'></a>\n",
    "$$\n",
    "h+d \\, (\\beta z^{-1})\\, d \\, (z) = 0 \\tag{13}\n",
    "$$\n",
    "\n",
    "Notice that if $ \\tilde z $ is a root of equation [(13)](#equation-oneseven), then so is $ \\beta \\tilde z^{-1} $.\n",
    "\n",
    "Thus, the roots of [(13)](#equation-oneseven) come in “$ \\beta $-reciprocal” pairs.\n",
    "\n",
    "Assume that the roots of [(13)](#equation-oneseven) are distinct.\n",
    "\n",
    "Let the roots be, in descending order according to their moduli, $ z_1, z_2, \\ldots, z_{2m} $.\n",
    "\n",
    "From the reciprocal pairs property and the assumption of distinct\n",
    "roots, it follows that $ \\vert z_j \\vert > \\sqrt \\beta\\ \\hbox{ for } j\\leq m \\hbox\n",
    "{ and } \\vert z_j \\vert < \\sqrt\\beta\\ \\hbox { for } j > m $.\n",
    "\n",
    "It also follows that $ z_{2m-j} = \\beta z^{-1}_{j+1}, j=0, 1, \\ldots, m-1 $.\n",
    "\n",
    "Therefore, the characteristic polynomial on the left side of [(13)](#equation-oneseven) can be expressed as\n",
    "\n",
    "\n",
    "<a id='equation-oneeight'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "h+d(\\beta z^{-1})d(z)\n",
    "&= z^{-m} z_0(z-z_1)\\cdots\n",
    "(z-z_m)(z-z_{m+1}) \\cdots (z-z_{2m}) \\cr\n",
    "&= z^{-m} z_0 (z-z_1)(z-z_2)\\cdots (z-z_m)(z-\\beta z_m^{-1})\n",
    "\\cdots  (z-\\beta z^{-1}_2)(z-\\beta z_1^{-1})\n",
    "\\end{aligned} \\tag{14}\n",
    "$$\n",
    "\n",
    "where $ z_0 $ is a constant.\n",
    "\n",
    "In [(14)](#equation-oneeight), we substitute $ (z-z_j) = -z_j (1- {1 \\over z_j}z) $ and\n",
    "$ (z-\\beta z_j^{-1}) = z(1 - {\\beta \\over z_j} z^{-1}) $ for $ j = 1, \\ldots, m $ to get\n",
    "\n",
    "$$\n",
    "h+d(\\beta z^{-1})d(z)\n",
    "= (-1)^m(z_0z_1\\cdots z_m)\n",
    "(1- {1\\over z_1} z) \\cdots (1-{1\\over z_m} z)(1- {1\\over z_1} \\beta z^{-1})\n",
    "\\cdots(1-{1\\over z_m} \\beta z^{-1})\n",
    "$$\n",
    "\n",
    "Now define $ c(z) = \\sum^m_{j=0} c_j \\, z^j $ as\n",
    "\n",
    "\n",
    "<a id='equation-onenine'></a>\n",
    "$$\n",
    "c\\,(z)=\\Bigl[(-1)^m z_0\\, z_1 \\cdots z_m\\Bigr]^{1/2} (1-{z\\over z_1}) \\,\n",
    "(1-{z\\over z_2}) \\cdots (1- {z\\over z_m}) \\tag{15}\n",
    "$$\n",
    "\n",
    "Notice that [(14)](#equation-oneeight) can be written\n",
    "\n",
    "\n",
    "<a id='equation-oneten'></a>\n",
    "$$\n",
    "h + d \\ (\\beta z^{-1})\\ d\\ (z) = c\\,(\\beta z^{-1})\\,c\\,(z) \\tag{16}\n",
    "$$\n",
    "\n",
    "It is useful to write [(15)](#equation-onenine) as\n",
    "\n",
    "\n",
    "<a id='equation-oneeleven'></a>\n",
    "$$\n",
    "c(z) = c_0(1-\\lambda_1\\, z) \\ldots (1-\\lambda_m z) \\tag{17}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "c_0\n",
    "= \\left[(-1)^m\\, z_0\\, z_1 \\cdots z_m\\right]^{1/2};\n",
    "\\quad \\lambda_j={1 \\over z_j},\\,\\  j=1, \\ldots, m\n",
    "$$\n",
    "\n",
    "Since $ \\vert z_j \\vert > \\sqrt \\beta \\hbox { for } j = 1, \\ldots, m $ it\n",
    "follows that $ \\vert \\lambda_j \\vert < 1/\\sqrt \\beta $ for $ j = 1,\n",
    "\\ldots, m $.\n",
    "\n",
    "Using [(17)](#equation-oneeleven), we can express the factorization [(16)](#equation-oneten) as\n",
    "\n",
    "$$\n",
    "h+d (\\beta z^{-1})d(z) = c^2_0 (1-\\lambda_1 z) \\cdots\n",
    "(1 - \\lambda_m z) (1-\\lambda_1 \\beta z^{-1})\n",
    "\\cdots (1 - \\lambda_m \\beta z^{-1})\n",
    "$$\n",
    "\n",
    "In sum, we have constructed a factorization [(16)](#equation-oneten) of the characteristic\n",
    "polynomial for the Euler equation in which the zeros of $ c(z) $\n",
    "exceed $ \\beta^{1/2} $ in modulus, and the zeros of\n",
    "$ c\\,(\\beta z^{-1}) $ are less than $ \\beta^{1/2} $ in modulus.\n",
    "\n",
    "Using [(16)](#equation-oneten), we now write the Euler equation as\n",
    "\n",
    "$$\n",
    "c(\\beta L^{-1})\\,c\\,(L)\\, y_t = a_t\n",
    "$$\n",
    "\n",
    "The unique solution of the Euler equation that satisfies condition [(12)](#equation-onesix)\n",
    "is\n",
    "\n",
    "\n",
    "<a id='equation-onethirteen'></a>\n",
    "$$\n",
    "c(L)\\,y_t = c\\,(\\beta L^{-1})^{-1}a_t \\tag{18}\n",
    "$$\n",
    "\n",
    "This can be established by using an argument paralleling that in\n",
    "chapter IX of [[Sar87]](https://python.quantecon.org/zreferences.html#sargent1987).\n",
    "\n",
    "To exhibit the solution in a form\n",
    "paralleling that of [[Sar87]](https://python.quantecon.org/zreferences.html#sargent1987), we use [(17)](#equation-oneeleven) to write\n",
    "[(18)](#equation-onethirteen) as\n",
    "\n",
    "\n",
    "<a id='equation-junk'></a>\n",
    "$$\n",
    "(1-\\lambda_1 L) \\cdots (1 - \\lambda_mL)y_t = {c^{-2}_0 a_t \\over (1-\\beta \\lambda_1 L^{-1}) \\cdots (1 - \\beta \\lambda_m L^{-1})} \\tag{19}\n",
    "$$\n",
    "\n",
    "Using [partial fractions](https://en.wikipedia.org/wiki/Partial_fraction_decomposition), we can write the characteristic polynomial on\n",
    "the right side of [(19)](#equation-junk) as\n",
    "\n",
    "$$\n",
    "\\sum^m_{j=1} {A_j \\over 1 - \\lambda_j \\, \\beta L^{-1}}\n",
    " \\quad \\text{where} \\quad\n",
    "A_j := {c^{-2}_0 \\over \\prod_{i \\not= j}(1-{\\lambda_i \\over \\lambda_j})}\n",
    "$$\n",
    "\n",
    "Then [(19)](#equation-junk) can be written\n",
    "\n",
    "$$\n",
    "(1-\\lambda_1 L) \\cdots (1-\\lambda_m L) y_t = \\sum^m_{j=1} \\, {A_j \\over 1 -\n",
    "\\lambda_j \\, \\beta L^{-1}} a_t\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "\n",
    "<a id='equation-onefifteen'></a>\n",
    "$$\n",
    "(1 - \\lambda_1 L) \\cdots (1 - \\lambda_m L) y_t = \\sum^m_{j=1}\\, A_j\n",
    "\\sum^\\infty_{k=0}\\, (\\lambda_j\\beta)^k\\, a_{t+k} \\tag{20}\n",
    "$$\n",
    "\n",
    "Equation [(20)](#equation-onefifteen) expresses the optimum sequence for $ y_t $ in terms\n",
    "of $ m $ lagged $ y $’s, and $ m $ weighted infinite\n",
    "geometric sums of future $ a_t $’s.\n",
    "\n",
    "Furthermore, [(20)](#equation-onefifteen) is the unique solution of the Euler equation that satisfies the initial conditions and condition [(12)](#equation-onesix).\n",
    "\n",
    "In effect, condition [(12)](#equation-onesix) compels us to\n",
    "solve the “unstable” roots of $ h+d (\\beta z^{-1})d(z) $ forward\n",
    "(see [[Sar87]](https://python.quantecon.org/zreferences.html#sargent1987)).\n",
    "\n",
    "The step of factoring the polynomial $ h+d (\\beta z^{-1})\\, d(z) $ into\n",
    "$ c\\, (\\beta z^{-1})c\\,(z) $, where the zeros of $ c\\,(z) $ all\n",
    "have modulus exceeding $ \\sqrt\\beta $, is central to solving the problem.\n",
    "\n",
    "We note two features of the solution [(20)](#equation-onefifteen)\n",
    "\n",
    "- Since $ \\vert \\lambda_j \\vert < 1/\\sqrt \\beta $ for all $ j $, it follows that $ (\\lambda_j \\ \\beta) < \\sqrt \\beta $.  \n",
    "- The assumption that $ \\{ a_t \\} $ is of exponential order less than $ 1 /\\sqrt \\beta $ is sufficient to guarantee that the geometric sums of future $ a_t $’s on the right side of [(20)](#equation-onefifteen) converge.  \n",
    "\n",
    "\n",
    "We immediately see that those sums will\n",
    "converge under the weaker condition that $ \\{ a_t\\} $ is of\n",
    "exponential order less than $ \\phi^{-1} $ where\n",
    "$ \\phi = \\max \\, \\{\\beta \\lambda_i, i=1,\\ldots,m\\} $.\n",
    "\n",
    "Note that with $ a_t $ identically zero, [(20)](#equation-onefifteen) implies that\n",
    "in general $ \\vert y_t \\vert $ eventually grows exponentially at a\n",
    "rate given by $ \\max_i \\vert \\lambda_i \\vert $.\n",
    "\n",
    "The condition\n",
    "$ \\max_i \\vert \\lambda_i \\vert <1 /\\sqrt \\beta $ guarantees that\n",
    "condition [(12)](#equation-onesix) is satisfied.\n",
    "\n",
    "In fact, $ \\max_i \\vert \\lambda_i\n",
    "\\vert < 1 /\\sqrt \\beta $ is a necessary condition for [(12)](#equation-onesix) to hold.\n",
    "\n",
    "Were [(12)](#equation-onesix) not satisfied, the objective function would diverge to $ - \\infty $, implying that the $ y_t $ path could not be optimal.\n",
    "\n",
    "For example, with $ a_t = 0 $, for all $ t $, it is easy to describe a naive (nonoptimal) policy for $ \\{y_t, t\\geq 0\\} $ that gives a finite value of [(17)](#equation-oneeleven).\n",
    "\n",
    "We can simply let $ y_t = 0 \\hbox { for } t\\geq 0 $.\n",
    "\n",
    "This policy involves at most $ m $ nonzero values of\n",
    "$ hy^2_t $ and $ [d(L)y_t]^2 $, and so yields a finite value of\n",
    "[(1)](#equation-oneone).\n",
    "\n",
    "Therefore it is easy to dominate a path that violates [(12)](#equation-onesix)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Undiscounted Problems\n",
    "\n",
    "It is worthwhile focusing on a special case of the LQ problems above:\n",
    "the undiscounted problem that emerges when $ \\beta = 1 $.\n",
    "\n",
    "In this case, the Euler equation is\n",
    "\n",
    "$$\n",
    "\\Bigl( h + d(L^{-1})d(L) \\Bigr)\\, y_t = a_t\n",
    "$$\n",
    "\n",
    "The factorization of the characteristic polynomial [(16)](#equation-oneten) becomes\n",
    "\n",
    "$$\n",
    "\\Bigl(h+d \\, (z^{-1})d(z)\\Bigr) = c\\,(z^{-1})\\, c\\,(z)\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "c\\,(z) &= c_0 (1 - \\lambda_1 z) \\ldots (1 - \\lambda_m z) \\cr\n",
    "c_0 &= \\Bigl[(-1)^m z_0 z_1 \\ldots z_m\\Bigr ] \\cr\n",
    "\\vert \\lambda_j \\vert &< 1 \\, \\hbox { for } \\, j = 1, \\ldots, m\\cr\n",
    "\\lambda_j &= \\frac{1}{z_j} \\hbox{ for } j=1,\\ldots, m\\cr\n",
    "z_0 &= \\hbox{ constant}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The solution of the problem becomes\n",
    "\n",
    "$$\n",
    "(1 - \\lambda_1 L) \\cdots (1 - \\lambda_m L) y_t = \\sum^m_{j=1} A_j\n",
    "\\sum^\\infty_{k=0} \\lambda^k_j a_{t+k}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transforming Discounted to Undiscounted Problem\n",
    "\n",
    "Discounted problems can always be converted into undiscounted problems via a simple transformation.\n",
    "\n",
    "Consider problem [(1)](#equation-oneone) with $ 0 < \\beta < 1 $.\n",
    "\n",
    "Define the transformed variables\n",
    "\n",
    "\n",
    "<a id='equation-onetwenty'></a>\n",
    "$$\n",
    "\\tilde a_t = \\beta^{t/2} a_t,\\ \\tilde y_t = \\beta^{t/2} y_t \\tag{21}\n",
    "$$\n",
    "\n",
    "Then notice that $ \\beta^t\\,[d\\, (L) y_t ]^2=[\\tilde d\\,(L)\\tilde y_t]^2 $ with\n",
    "$ \\tilde d \\,(L)=\\sum^m_{j=0} \\tilde d_j\\, L^j $ and $ \\tilde d_j = \\beta^{j/2}\n",
    "d_j $.\n",
    "\n",
    "Then the original criterion function [(1)](#equation-oneone) is equivalent to\n",
    "\n",
    "\n",
    "<a id='equation-oneoneprime'></a>\n",
    "$$\n",
    "\\lim_{N \\rightarrow \\infty}\n",
    "\\sum^N_{t=0}\n",
    "\\{\\tilde a_t\\, \\tilde y_t - {1 \\over 2} h\\,\\tilde y^2_t - {1\\over 2}\n",
    "[ \\tilde d\\,(L)\\, \\tilde y_t]^2 \\} \\tag{22}\n",
    "$$\n",
    "\n",
    "which is to be maximized over sequences $ \\{\\tilde y_t,\\  t=0, \\ldots\\} $ subject to\n",
    "$ \\tilde y_{-1}, \\cdots, \\tilde y_{-m} $ given and $ \\{\\tilde a_t,\\  t=1, \\ldots\\} $ a known bounded sequence.\n",
    "\n",
    "The Euler equation for this problem is $ [h+\\tilde d \\,(L^{-1}) \\, \\tilde d\\, (L) ]\\, \\tilde y_t = \\tilde a_t $.\n",
    "\n",
    "The solution is\n",
    "\n",
    "$$\n",
    "(1 - \\tilde \\lambda_1 L) \\cdots (1 - \\tilde \\lambda_m L)\\,\\tilde y_t =\n",
    "\\sum^m_{j=1} \\tilde A_j \\sum^\\infty_{k=0} \\tilde \\lambda^k_j \\, \\tilde a_{t+k}\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "\n",
    "<a id='equation-onetwentyone'></a>\n",
    "$$\n",
    "\\tilde y_t = \\tilde f_1 \\, \\tilde y_{t-1} + \\cdots + \\tilde f_m\\,\n",
    "\\tilde y_{t-m} + \\sum^m_{j=1} \\tilde A_j \\sum^\\infty_{k=0} \\tilde \\lambda^k_j\n",
    "\\, \\tilde a_{t+k}, \\tag{23}\n",
    "$$\n",
    "\n",
    "where $ \\tilde c \\,(z^{-1}) \\tilde c\\,(z) = h + \\tilde d\\,(z^{-1}) \\tilde d \\,(z) $, and where\n",
    "\n",
    "$$\n",
    "\\bigl[(-1)^m\\, \\tilde z_0 \\tilde z_1 \\ldots \\tilde z_m \\bigr]^{1/2}\n",
    "(1 - \\tilde \\lambda_1\\, z) \\ldots (1 - \\tilde \\lambda_m\\, z) = \\tilde c\\,(z),\n",
    "\\hbox { where } \\ \\vert \\tilde \\lambda_j \\vert < 1\n",
    "$$\n",
    "\n",
    "We leave it to the reader to show that [(23)](#equation-onetwentyone) implies the equivalent form of the solution\n",
    "\n",
    "$$\n",
    "y_t = f_1\\, y_{t-1} + \\cdots + f_m\\, y_{t-m} + \\sum^m_{j=1} A_j\n",
    "\\sum^\\infty_{k=0} \\, (\\lambda_j\\, \\beta)^k \\, a_{t+k}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "\n",
    "<a id='equation-onetwentythree'></a>\n",
    "$$\n",
    "f_j = \\tilde f_j\\, \\beta^{-j/2},\\  A_j = \\tilde A_j,\\  \\lambda_j = \\tilde\n",
    "\\lambda_j \\, \\beta^{-1/2} \\tag{24}\n",
    "$$\n",
    "\n",
    "The transformations [(21)](#equation-onetwenty) and the inverse formulas\n",
    "[(24)](#equation-onetwentythree) allow us to solve a discounted problem by first\n",
    "solving a related undiscounted problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "\n",
    "Here’s the code that computes solutions to the LQ problem using the methods described\n",
    "above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as spst\n",
    "import scipy.linalg as la\n",
    "\n",
    "class LQFilter:\n",
    "\n",
    "    def __init__(self, d, h, y_m, r=None, h_eps=None, β=None):\n",
    "        \"\"\"\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "            d : list or numpy.array (1-D or a 2-D column vector)\n",
    "                    The order of the coefficients: [d_0, d_1, ..., d_m]\n",
    "            h : scalar\n",
    "                    Parameter of the objective function (corresponding to the\n",
    "                    quadratic term)\n",
    "            y_m : list or numpy.array (1-D or a 2-D column vector)\n",
    "                    Initial conditions for y\n",
    "            r : list or numpy.array (1-D or a 2-D column vector)\n",
    "                    The order of the coefficients: [r_0, r_1, ..., r_k]\n",
    "                    (optional, if not defined -> deterministic problem)\n",
    "            β : scalar\n",
    "                    Discount factor (optional, default value is one)\n",
    "        \"\"\"\n",
    "\n",
    "        self.h = h\n",
    "        self.d = np.asarray(d)\n",
    "        self.m = self.d.shape[0] - 1\n",
    "\n",
    "        self.y_m = np.asarray(y_m)\n",
    "\n",
    "        if self.m == self.y_m.shape[0]:\n",
    "            self.y_m = self.y_m.reshape(self.m, 1)\n",
    "        else:\n",
    "            raise ValueError(\"y_m must be of length m = {self.m:d}\")\n",
    "\n",
    "        #---------------------------------------------\n",
    "        # Define the coefficients of ϕ upfront\n",
    "        #---------------------------------------------\n",
    "        ϕ = np.zeros(2 * self.m + 1)\n",
    "        for i in range(- self.m, self.m + 1):\n",
    "            ϕ[self.m - i] = np.sum(np.diag(self.d.reshape(self.m + 1, 1) \\\n",
    "                                           @ self.d.reshape(1, self.m + 1),\n",
    "                                           k=-i\n",
    "                                           )\n",
    "                                    )\n",
    "        ϕ[self.m] = ϕ[self.m] + self.h\n",
    "        self.ϕ = ϕ\n",
    "\n",
    "        #-----------------------------------------------------\n",
    "        # If r is given calculate the vector ϕ_r\n",
    "        #-----------------------------------------------------\n",
    "        if r is None:\n",
    "            pass\n",
    "        else:\n",
    "            self.r = np.asarray(r)\n",
    "            self.k = self.r.shape[0] - 1\n",
    "            ϕ_r = np.zeros(2 * self.k + 1)\n",
    "            for i in range(- self.k, self.k + 1):\n",
    "                ϕ_r[self.k - i] = np.sum(np.diag(self.r.reshape(self.k + 1, 1) \\\n",
    "                                                 @ self.r.reshape(1, self.k + 1),\n",
    "                                                 k=-i\n",
    "                                                 )\n",
    "                                        )\n",
    "            if h_eps is None:\n",
    "                self.ϕ_r = ϕ_r\n",
    "            else:\n",
    "                ϕ_r[self.k] = ϕ_r[self.k] + h_eps\n",
    "                self.ϕ_r = ϕ_r\n",
    "\n",
    "        #-----------------------------------------------------\n",
    "        # If β is given, define the transformed variables\n",
    "        #-----------------------------------------------------\n",
    "        if β is None:\n",
    "            self.β = 1\n",
    "        else:\n",
    "            self.β = β\n",
    "            self.d = self.β**(np.arange(self.m + 1)/2) * self.d\n",
    "            self.y_m = self.y_m * (self.β**(- np.arange(1, self.m + 1)/2)) \\\n",
    "                                   .reshape(self.m, 1)\n",
    "\n",
    "    def construct_W_and_Wm(self, N):\n",
    "        \"\"\"\n",
    "        This constructs the matrices W and W_m for a given number of periods N\n",
    "        \"\"\"\n",
    "\n",
    "        m = self.m\n",
    "        d = self.d\n",
    "\n",
    "        W = np.zeros((N + 1, N + 1))\n",
    "        W_m = np.zeros((N + 1, m))\n",
    "\n",
    "        #---------------------------------------\n",
    "        # Terminal conditions\n",
    "        #---------------------------------------\n",
    "\n",
    "        D_m1 = np.zeros((m + 1, m + 1))\n",
    "        M = np.zeros((m + 1, m))\n",
    "\n",
    "        # (1) Constuct the D_{m+1} matrix using the formula\n",
    "\n",
    "        for j in range(m + 1):\n",
    "            for k in range(j, m + 1):\n",
    "                D_m1[j, k] = d[:j + 1] @ d[k - j: k + 1]\n",
    "\n",
    "        # Make the matrix symmetric\n",
    "        D_m1 = D_m1 + D_m1.T - np.diag(np.diag(D_m1))\n",
    "\n",
    "        # (2) Construct the M matrix using the entries of D_m1\n",
    "\n",
    "        for j in range(m):\n",
    "            for i in range(j + 1, m + 1):\n",
    "                M[i, j] = D_m1[i - j - 1, m]\n",
    "\n",
    "        #----------------------------------------------\n",
    "        # Euler equations for t = 0, 1, ..., N-(m+1)\n",
    "        #----------------------------------------------\n",
    "        ϕ = self.ϕ\n",
    "\n",
    "        W[:(m + 1), :(m + 1)] = D_m1 + self.h * np.eye(m + 1)\n",
    "        W[:(m + 1), (m + 1):(2 * m + 1)] = M\n",
    "\n",
    "        for i, row in enumerate(np.arange(m + 1, N + 1 - m)):\n",
    "            W[row, (i + 1):(2 * m + 2 + i)] = ϕ\n",
    "\n",
    "        for i in range(1, m + 1):\n",
    "            W[N - m + i, -(2 * m + 1 - i):] = ϕ[:-i]\n",
    "\n",
    "        for i in range(m):\n",
    "            W_m[N - i, :(m - i)] = ϕ[(m + 1 + i):]\n",
    "\n",
    "        return W, W_m\n",
    "\n",
    "    def roots_of_characteristic(self):\n",
    "        \"\"\"\n",
    "        This function calculates z_0 and the 2m roots of the characteristic\n",
    "        equation associated with the Euler equation (1.7)\n",
    "\n",
    "        Note:\n",
    "        ------\n",
    "        numpy.poly1d(roots, True) defines a polynomial using its roots that can\n",
    "        be evaluated at any point. If x_1, x_2, ... , x_m are the roots then\n",
    "            p(x) = (x - x_1)(x - x_2)...(x - x_m)\n",
    "        \"\"\"\n",
    "        m = self.m\n",
    "        ϕ = self.ϕ\n",
    "\n",
    "        # Calculate the roots of the 2m-polynomial\n",
    "        roots = np.roots(ϕ)\n",
    "        # Sort the roots according to their length (in descending order)\n",
    "        roots_sorted = roots[np.argsort(abs(roots))[::-1]]\n",
    "\n",
    "        z_0 = ϕ.sum() / np.poly1d(roots, True)(1)\n",
    "        z_1_to_m = roots_sorted[:m]     # We need only those outside the unit circle\n",
    "\n",
    "        λ = 1 / z_1_to_m\n",
    "\n",
    "        return z_1_to_m, z_0, λ\n",
    "\n",
    "    def coeffs_of_c(self):\n",
    "        '''\n",
    "        This function computes the coefficients {c_j, j = 0, 1, ..., m} for\n",
    "                c(z) = sum_{j = 0}^{m} c_j z^j\n",
    "\n",
    "        Based on the expression (1.9). The order is\n",
    "            c_coeffs = [c_0, c_1, ..., c_{m-1}, c_m]\n",
    "        '''\n",
    "        z_1_to_m, z_0 = self.roots_of_characteristic()[:2]\n",
    "\n",
    "        c_0 = (z_0 * np.prod(z_1_to_m).real * (- 1)**self.m)**(.5)\n",
    "        c_coeffs = np.poly1d(z_1_to_m, True).c * z_0 / c_0\n",
    "\n",
    "        return c_coeffs[::-1]\n",
    "\n",
    "    def solution(self):\n",
    "        \"\"\"\n",
    "        This function calculates {λ_j, j=1,...,m} and {A_j, j=1,...,m}\n",
    "        of the expression (1.15)\n",
    "        \"\"\"\n",
    "        λ = self.roots_of_characteristic()[2]\n",
    "        c_0 = self.coeffs_of_c()[-1]\n",
    "\n",
    "        A = np.zeros(self.m, dtype=complex)\n",
    "        for j in range(self.m):\n",
    "            denom = 1 - λ/λ[j]\n",
    "            A[j] = c_0**(-2) / np.prod(denom[np.arange(self.m) != j])\n",
    "\n",
    "        return λ, A\n",
    "\n",
    "    def construct_V(self, N):\n",
    "        '''\n",
    "        This function constructs the covariance matrix for x^N (see section 6)\n",
    "        for a given period N\n",
    "        '''\n",
    "        V = np.zeros((N, N))\n",
    "        ϕ_r = self.ϕ_r\n",
    "\n",
    "        for i in range(N):\n",
    "            for j in range(N):\n",
    "                if abs(i-j) <= self.k:\n",
    "                    V[i, j] = ϕ_r[self.k + abs(i-j)]\n",
    "\n",
    "        return V\n",
    "\n",
    "    def simulate_a(self, N):\n",
    "        \"\"\"\n",
    "        Assuming that the u's are normal, this method draws a random path\n",
    "        for x^N\n",
    "        \"\"\"\n",
    "        V = self.construct_V(N + 1)\n",
    "        d = spst.multivariate_normal(np.zeros(N + 1), V)\n",
    "\n",
    "        return d.rvs()\n",
    "\n",
    "    def predict(self, a_hist, t):\n",
    "        \"\"\"\n",
    "        This function implements the prediction formula discussed in section 6 (1.59)\n",
    "        It takes a realization for a^N, and the period in which the prediction is \n",
    "        formed\n",
    "\n",
    "        Output:  E[abar | a_t, a_{t-1}, ..., a_1, a_0]\n",
    "        \"\"\"\n",
    "\n",
    "        N = np.asarray(a_hist).shape[0] - 1\n",
    "        a_hist = np.asarray(a_hist).reshape(N + 1, 1)\n",
    "        V = self.construct_V(N + 1)\n",
    "\n",
    "        aux_matrix = np.zeros((N + 1, N + 1))\n",
    "        aux_matrix[:(t + 1), :(t + 1)] = np.eye(t + 1)\n",
    "        L = la.cholesky(V).T\n",
    "        Ea_hist = la.inv(L) @ aux_matrix @ L @ a_hist\n",
    "\n",
    "        return Ea_hist\n",
    "\n",
    "    def optimal_y(self, a_hist, t=None):\n",
    "        \"\"\"\n",
    "        - if t is NOT given it takes a_hist (list or numpy.array) as a\n",
    "          deterministic a_t\n",
    "        - if t is given, it solves the combined control prediction problem \n",
    "          (section 7)(by default, t == None -> deterministic)\n",
    "\n",
    "        for a given sequence of a_t (either deterministic or a particular\n",
    "        realization), it calculates the optimal y_t sequence using the method\n",
    "        of the lecture\n",
    "\n",
    "        Note:\n",
    "        ------\n",
    "        scipy.linalg.lu normalizes L, U so that L has unit diagonal elements\n",
    "        To make things consistent with the lecture, we need an auxiliary\n",
    "        diagonal matrix D which renormalizes L and U\n",
    "        \"\"\"\n",
    "\n",
    "        N = np.asarray(a_hist).shape[0] - 1\n",
    "        W, W_m = self.construct_W_and_Wm(N)\n",
    "\n",
    "        L, U = la.lu(W, permute_l=True)\n",
    "        D = np.diag(1 / np.diag(U))\n",
    "        U = D @ U\n",
    "        L = L @ np.diag(1 / np.diag(D))\n",
    "\n",
    "        J = np.fliplr(np.eye(N + 1))\n",
    "\n",
    "        if t is None:   # If the problem is deterministic\n",
    "\n",
    "            a_hist = J @ np.asarray(a_hist).reshape(N + 1, 1)\n",
    "\n",
    "            #--------------------------------------------\n",
    "            # Transform the 'a' sequence if β is given\n",
    "            #--------------------------------------------\n",
    "            if self.β != 1:\n",
    "                a_hist =  a_hist * (self.β**(np.arange(N + 1) / 2))[::-1] \\\n",
    "                                    .reshape(N + 1, 1)\n",
    "\n",
    "            a_bar = a_hist - W_m @ self.y_m            # a_bar from the lecture\n",
    "            Uy = np.linalg.solve(L, a_bar)             # U @ y_bar = L^{-1}\n",
    "            y_bar = np.linalg.solve(U, Uy)             # y_bar = U^{-1}L^{-1}\n",
    "\n",
    "            # Reverse the order of y_bar with the matrix J\n",
    "            J = np.fliplr(np.eye(N + self.m + 1))\n",
    "            # y_hist : concatenated y_m and y_bar\n",
    "            y_hist = J @ np.vstack([y_bar, self.y_m])\n",
    "\n",
    "            #--------------------------------------------\n",
    "            # Transform the optimal sequence back if β is given\n",
    "            #--------------------------------------------\n",
    "            if self.β != 1:\n",
    "                y_hist = y_hist * (self.β**(- np.arange(-self.m, N + 1)/2)) \\\n",
    "                                    .reshape(N + 1 + self.m, 1)\n",
    "\n",
    "            return y_hist, L, U, y_bar\n",
    "\n",
    "        else:           # If the problem is stochastic and we look at it\n",
    "\n",
    "            Ea_hist = self.predict(a_hist, t).reshape(N + 1, 1)\n",
    "            Ea_hist = J @ Ea_hist\n",
    "\n",
    "            a_bar = Ea_hist - W_m @ self.y_m           # a_bar from the lecture\n",
    "            Uy = np.linalg.solve(L, a_bar)             # U @ y_bar = L^{-1}\n",
    "            y_bar = np.linalg.solve(U, Uy)             # y_bar = U^{-1}L^{-1}\n",
    "\n",
    "            # Reverse the order of y_bar with the matrix J\n",
    "            J = np.fliplr(np.eye(N + self.m + 1))\n",
    "            # y_hist : concatenated y_m and y_bar\n",
    "            y_hist = J @ np.vstack([y_bar, self.y_m])\n",
    "\n",
    "            return y_hist, L, U, y_bar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "In this application, we’ll have one lag, with\n",
    "\n",
    "$$\n",
    "d(L) y_t =  \\gamma(I - L) y_t = \\gamma (y_t - y_{t-1})\n",
    "$$\n",
    "\n",
    "Suppose for the moment that $ \\gamma = 0 $.\n",
    "\n",
    "Then the intertemporal component of the LQ problem disappears, and the agent\n",
    "simply wants to maximize $ a_t y_t - hy^2_t / 2 $ in each period.\n",
    "\n",
    "This means that the agent chooses $ y_t = a_t / h $.\n",
    "\n",
    "In the following we’ll set $ h = 1 $, so that the agent just wants to\n",
    "track the $ \\{a_t\\} $ process.\n",
    "\n",
    "However, as we increase $ \\gamma $, the agent gives greater weight to a smooth time path.\n",
    "\n",
    "Hence $ \\{y_t\\} $ evolves as a smoothed version of $ \\{a_t\\} $.\n",
    "\n",
    "The $ \\{a_t\\} $ sequence we’ll choose as a stationary cyclic process plus some white noise.\n",
    "\n",
    "Here’s some code that generates a plot when $ \\gamma = 0.8 $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Set seed and generate a_t sequence\n",
    "np.random.seed(123)\n",
    "n = 100\n",
    "a_seq = np.sin(np.linspace(0, 5 * np.pi, n)) + 2 + 0.1 * np.random.randn(n)\n",
    "\n",
    "def plot_simulation(γ=0.8, m=1, h=1, y_m=2):\n",
    "\n",
    "    d = γ * np.asarray([1, -1])\n",
    "    y_m = np.asarray(y_m).reshape(m, 1)\n",
    "\n",
    "    testlq = LQFilter(d, h, y_m)\n",
    "    y_hist, L, U, y = testlq.optimal_y(a_seq)\n",
    "    y = y[::-1]  # Reverse y\n",
    "\n",
    "    # Plot simulation results\n",
    "\n",
    "    fig, ax = plt.subplots(figsize=(10, 6))\n",
    "    p_args = {'lw' : 2, 'alpha' : 0.6}\n",
    "    time = range(len(y))\n",
    "    ax.plot(time, a_seq / h, 'k-o', ms=4, lw=2, alpha=0.6, label='$a_t$')\n",
    "    ax.plot(time, y, 'b-o', ms=4, lw=2, alpha=0.6, label='$y_t$')\n",
    "    ax.set(title=rf'Dynamics with $\\gamma = {γ}$',\n",
    "           xlabel='Time',\n",
    "           xlim=(0, max(time))\n",
    "          )\n",
    "    ax.legend()\n",
    "    ax.grid()\n",
    "    plt.show()\n",
    "\n",
    "plot_simulation()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here’s what happens when we change $ \\gamma $ to 5.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "plot_simulation(γ=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here’s $ \\gamma = 10 $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "plot_simulation(γ=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Consider solving a discounted version $ (\\beta < 1) $ of problem\n",
    "[(1)](#equation-oneone), as follows.\n",
    "\n",
    "Convert [(1)](#equation-oneone) to the undiscounted problem [(22)](#equation-oneoneprime).\n",
    "\n",
    "Let the solution of [(22)](#equation-oneoneprime) in feedback form be\n",
    "\n",
    "$$\n",
    "(1 - \\tilde \\lambda_1 L)\\, \\cdots\\, (1 - \\tilde \\lambda_m L) \\tilde y_t =\n",
    "\\sum^m_{j=1} \\tilde A_j \\sum^\\infty_{k=0} \\tilde \\lambda^k_j \\tilde a_{t+k}\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "\n",
    "<a id='equation-estar'></a>\n",
    "$$\n",
    "\\tilde y_t = \\tilde f_1 \\tilde y_{t-1} + \\cdots + \\tilde f_m \\tilde y_{t-m} +\n",
    "\\sum^m_{j=1} \\tilde A_j \\sum^\\infty_{k=0} \\tilde \\lambda^k_j \\tilde a_{t+k} \\tag{25}\n",
    "$$\n",
    "\n",
    "Here\n",
    "\n",
    "- $ h + \\tilde d (z^{-1}) \\tilde d (z) = \\tilde c (z^{-1}) \\tilde c (z) $  \n",
    "- $ \\tilde c (z) = [(-1)^m \\tilde z_0 \\tilde z_1 \\cdots \\tilde z_m ]^{1/2} (1 - \\tilde \\lambda_1 z) \\cdots (1 - \\tilde \\lambda_m z) $  \n",
    "\n",
    "\n",
    "where the $ \\tilde z_j $ are the zeros of $ h +\\tilde d (z^{-1})\\, \\tilde d(z) $.\n",
    "\n",
    "Prove that [(25)](#equation-estar) implies that the solution for $ y_t $ in feedback form is\n",
    "\n",
    "$$\n",
    "y_t = f_1 y_{t-1} + \\ldots + f_m y_{t-m} + \\sum^m_{j=1} A_j\n",
    "\\sum^\\infty_{k=0} \\beta^k \\lambda^k_j a_{t+k}\n",
    "$$\n",
    "\n",
    "where $ f_j = \\tilde f_j \\beta^{-j/2}, A_j = \\tilde A_j $, and $ \\lambda_j = \\tilde \\lambda_j \\beta^{-1/2} $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "\n",
    "Solve the optimal control problem, maximize\n",
    "\n",
    "$$\n",
    "\\sum^2_{t=0}\\ \\Bigl\\{a_t y_t - {1 \\over 2} [(1 - 2 L) y_t]^2\\Bigr\\}\n",
    "$$\n",
    "\n",
    "subject to $ y_{-1} $ given, and $ \\{ a_t\\} $ a known bounded sequence.\n",
    "\n",
    "Express the solution in the “feedback form” [(20)](#equation-onefifteen), giving numerical values for the coefficients.\n",
    "\n",
    "Make sure that the boundary conditions [(5)](#equation-onefive) are satisfied.\n",
    "\n",
    "(Note: this problem differs from the problem in the text in one important way: instead of $ h > 0 $ in [(1)](#equation-oneone), $ h = 0 $. This has an important influence on the solution.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3\n",
    "\n",
    "Solve the infinite time-optimal control problem to maximize\n",
    "\n",
    "$$\n",
    "\\lim_{N \\rightarrow \\infty}\n",
    "\\sum^N_{t=0}\\, -\\, {1 \\over 2} [(1 -2 L) y_t]^2,\n",
    "$$\n",
    "\n",
    "subject to $ y_{-1} $ given. Prove that the solution is\n",
    "\n",
    "$$\n",
    "y_t = 2y_{t-1} = 2^{t+1} y_{-1} \\qquad t > 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4\n",
    "\n",
    "Solve the infinite time problem, to maximize\n",
    "\n",
    "$$\n",
    "\\lim_{N \\rightarrow \\infty}\\ \\sum^N_{t=0}\\ (.0000001)\\, y^2_t - {1 \\over 2}\n",
    "[(1 - 2 L) y_t]^2\n",
    "$$\n",
    "\n",
    "subject to $ y_{-1} $ given. Prove that the solution $ y_t = 2y_{t-1} $  violates condition [(12)](#equation-onesix), and so\n",
    "is not optimal.\n",
    "\n",
    "Prove that the optimal solution is approximately $ y_t = .5 y_{t-1} $."
   ]
  }
 ],
 "metadata": {
  "date": 1584334744.6397896,
  "filename": "lu_tricks.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Classical Control with Linear Algebra"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}