lagrangian_lqdp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b4c34f88",
   "metadata": {},
   "source": [
    "# Lagrangian for LQ Control"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "865acb46",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "!pip install quantecon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "548ade1e",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from quantecon import LQ\n",
    "from scipy.linalg import schur"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40618c64",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "This is a sequel to this lecture [linear quadratic dynamic programming](https://python.quantecon.org/lqcontrol.html)\n",
    "\n",
    "It can also be regarded as presenting **invariant subspace** techniques that extend ones that\n",
    "we encountered earlier in this lecture [stability in linear rational expectations models](https://python.quantecon.org/re_with_feedback.html)\n",
    "\n",
    "We present a Lagrangian formulation of an infinite horizon linear quadratic undiscounted dynamic programming problem.\n",
    "\n",
    "Such a problem is  also sometimes called an  optimal linear regulator problem.\n",
    "\n",
    "A Lagrangian formulation\n",
    "\n",
    "- carries insights about connections between stability and optimality  \n",
    "- is the basis for fast algorithms for solving Riccati equations  \n",
    "- opens the way to constructing solutions of dynamic systems that don’t come directly from an  intertemporal optimization problem  \n",
    "\n",
    "\n",
    "A key tool in this lecture is the concept of an $ n \\times n $ **symplectic** matrix.\n",
    "\n",
    "A symplectic matrix has eigenvalues that occur in **reciprocal pairs**, meaning that if $ \\lambda_i \\in (-1,1) $ is an eigenvalue, then so is $ \\lambda_i^{-1} $.\n",
    "\n",
    "This reciprocal pairs property of the eigenvalues of a matrix  is a tell-tale sign that the matrix describes the joint dynamics of a system of\n",
    "equations describing the **states** and **costates** that constitute   first-order necessary  conditions for solving an undiscounted linear-quadratic  infinite-horizon optimization problem.\n",
    "\n",
    "The symplectic matrix that will interest us describes the first-order dynamics of **state** and **co-state** vectors of an optimally controlled system.\n",
    "\n",
    "In focusing on eigenvalues and eigenvectors of this matrix, we capitalize on an analysis of\n",
    "**invariant subspaces.**\n",
    "\n",
    "These invariant subspace formulations of LQ dynamic programming problems provide a bridge between recursive\n",
    "(i.e., dynamic programming) formulations and classical formulations of linear control and linear filtering problems that make use of related matrix decompositions (see for example [this lecture](https://python-advanced.quantecon.org/lu_tricks.html) and [this lecture](https://python-advanced.quantecon.org/classical_filtering.html)).\n",
    "\n",
    "While most of this lecture focuses on undiscounted problems, later sections describe handy ways of transforming discounted problems to undiscounted ones.\n",
    "\n",
    "The techniques in this lecture will prove useful when we study Stackelberg and Ramsey problem in\n",
    "[this lecture](https://python-advanced.quantecon.org/dyn_stack.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc32f8bb",
   "metadata": {},
   "source": [
    "## Undiscounted LQ DP Problem\n",
    "\n",
    "The problem is to choose a sequence of controls  $ \\{u_t\\}_{t=0}^\\infty $ to maximize the criterion\n",
    "\n",
    "$$\n",
    "- \\sum_{t=0}^\\infty \\{x'_t Rx_t+u'_tQu_t\\}\n",
    "$$\n",
    "\n",
    "subject to $ x_{t+1}=Ax_t+Bu_t $, where  $ x_0 $ is a  given initial state vector.\n",
    "\n",
    "Here $ x_t $ is an $ (n\\times 1) $ vector of state variables, $ u_t $ is a $ (k\\times 1) $\n",
    "vector of controls, $ R $ is a positive semidefinite symmetric matrix,\n",
    "$ Q $ is a positive definite symmetric matrix, $ A $ is an $ (n\\times n) $\n",
    "matrix, and $ B $ is an $ (n\\times k) $ matrix.\n",
    "\n",
    "The optimal\n",
    "value function  turns out to be  quadratic, $ V(x)= - x'Px $, where $ P $ is a positive\n",
    "semidefinite symmetric matrix.\n",
    "\n",
    "Using the transition law to eliminate next period’s state, the Bellman\n",
    "equation becomes\n",
    "\n",
    "\n",
    "<a id='equation-bellman0'></a>\n",
    "$$\n",
    "- x'Px=\\max_u \\{- x' Rx-u'Qu-(Ax+Bu)' P(Ax+Bu)\\} \\tag{57.1}\n",
    "$$\n",
    "\n",
    "The first-order necessary conditions for the maximum problem on the\n",
    "right side of equation [(57.1)](#equation-bellman0) are\n",
    "\n",
    ">**Note**\n",
    ">\n",
    ">We use the following rules for differentiating quadratic and bilinear matrix forms:\n",
    "$ {\\partial x' A x \\over \\partial x} = (A + A') x; {\\partial y' B z \\over \\partial y} = B z, {\\partial\n",
    "y' B z \\over \\partial z} = B' y $.\n",
    "\n",
    "$$\n",
    "(Q+B'PB)u=-B'PAx,\n",
    "$$\n",
    "\n",
    "which implies that an optimal decision rule for $ u $ is\n",
    "\n",
    "$$\n",
    "u=-(Q+B'PB)^{-1} B'PAx\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "$$\n",
    "u=-Fx,\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "F=(Q+B'PB)^{-1}B'PA.\n",
    "$$\n",
    "\n",
    "Substituting $ u = - (Q+B'PB)^{-1}B'PAx $  into\n",
    "the right side of equation [(57.1)](#equation-bellman0) and rearranging gives\n",
    "\n",
    "\n",
    "<a id='equation-riccati'></a>\n",
    "$$\n",
    "P=R+A'PA-A'PB(Q+B'PB)^{-1} B'PA. \\tag{57.2}\n",
    "$$\n",
    "\n",
    "Equation [(57.2)](#equation-riccati) is called an **algebraic matrix Riccati** equation.\n",
    "\n",
    "There are multiple solutions of equation [(57.2)](#equation-riccati).\n",
    "\n",
    "But only one of them is positive definite.\n",
    "\n",
    "The positive define solution is associated with the maximum of our problem.\n",
    "\n",
    "It expresses the matrix $ P $ as an implicit function of the matrices\n",
    "$ R,Q,A,B $.\n",
    "\n",
    "Notice that the **gradient of the value function** is\n",
    "\n",
    "\n",
    "<a id='equation-eqn-valgrad'></a>\n",
    "$$\n",
    "\\frac{\\partial V(x)}{\\partial x} = - 2 P x \\tag{57.3}\n",
    "$$\n",
    "\n",
    "We shall use fact [(57.3)](#equation-eqn-valgrad) later."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddf3723f",
   "metadata": {},
   "source": [
    "## Lagrangian\n",
    "\n",
    "For the undiscounted optimal linear regulator problem, form the Lagrangian\n",
    "\n",
    "\n",
    "<a id='equation-lag-lqdp-eq1'></a>\n",
    "$$\n",
    "{\\cal L} = - \\sum^\\infty_{t=0} \\biggl\\{ x^\\prime_t R x_t + u_t^\\prime Q u_t +\n",
    "                                 2 \\mu^\\prime_{t+1} [A x_t + B u_t - x_{t+1}]\\biggr\\} \\tag{57.4}\n",
    "$$\n",
    "\n",
    "where $ 2 \\mu_{t+1} $ is a vector of Lagrange multipliers on the time $ t $ transition law $ x_{t+1} = A x_t + B u_t $.\n",
    "\n",
    "(We put the $ 2 $ in front of $ \\mu_{t+1} $ to make things match up nicely with equation [(57.3)](#equation-eqn-valgrad).)\n",
    "\n",
    "First-order conditions for maximization with respect to $ \\{u_t,x_{t+1}\\}_{t=0}^\\infty $ are\n",
    "\n",
    "\n",
    "<a id='equation-lag-lqdp-eq2'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "2 Q u_t &+ 2B^\\prime \\mu_{t+1} = 0 \\ ,\\ t \\geq 0 \\cr \\mu_t &= R x_t + A^\\prime \\mu_{t+1}\\ ,\\ t\\geq 1.\\cr\n",
    "\\end{aligned} \\tag{57.5}\n",
    "$$\n",
    "\n",
    "Define $ \\mu_0 $ to be  a vector of shadow prices of $ x_0 $ and apply an envelope condition to [(57.4)](#equation-lag-lqdp-eq1)\n",
    "to deduce that\n",
    "\n",
    "$$\n",
    "\\mu_0 = R x_0 + A' \\mu_1,\n",
    "$$\n",
    "\n",
    "which is a time $ t=0 $ counterpart to the second equation of system [(57.5)](#equation-lag-lqdp-eq2).\n",
    "\n",
    "An important fact is  that\n",
    "\n",
    "\n",
    "<a id='equation-eqn-mupx'></a>\n",
    "$$\n",
    "\\mu_{t+1} = P x_{t+1} \\tag{57.6}\n",
    "$$\n",
    "\n",
    "where $ P $ is a positive define  matrix that solves  the algebraic Riccati equation [(57.2)](#equation-riccati).\n",
    "\n",
    "Thus, from equations [(57.3)](#equation-eqn-valgrad) and  [(57.6)](#equation-eqn-mupx),  $ - 2 \\mu_{t} $ is\n",
    "the gradient of the value function with respect to $ x_t $.\n",
    "\n",
    "The Lagrange multiplier vector $ \\mu_{t} $ is often called the **costate** vector that\n",
    "corresponds to the **state** vector $ x_t $.\n",
    "\n",
    "It is useful to proceed with the following steps:\n",
    "\n",
    "- solve the first equation of [(57.5)](#equation-lag-lqdp-eq2)  for $ u_t $ in terms of $ \\mu_{t+1} $.  \n",
    "- substitute the result into the law of motion $ x_{t+1} = A x_t + B u_t $.  \n",
    "- arrange the resulting equation and the second equation of [(57.5)](#equation-lag-lqdp-eq2)  into the form  \n",
    "\n",
    "\n",
    "\n",
    "<a id='equation-eq-systosolve'></a>\n",
    "$$\n",
    "L\\ \\begin{pmatrix}x_{t+1}\\cr \\mu_{t+1}\\cr\\end{pmatrix}\\ = \\ N\\ \\begin{pmatrix}x_t\\cr \\mu_t\\cr\\end{pmatrix}\\\n",
    ",\\ t \\geq 0, \\tag{57.7}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "L = \\ \\begin{pmatrix}I & BQ^{-1} B^\\prime \\cr 0 & A^\\prime\\cr\\end{pmatrix}, \\quad N = \\\n",
    "\\begin{pmatrix}A & 0\\cr -R & I\\cr\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "When $ L $ is of full rank (i.e., when $ A $ is of full rank), we can write\n",
    "system [(57.7)](#equation-eq-systosolve) as\n",
    "\n",
    "\n",
    "<a id='equation-eq4orig'></a>\n",
    "$$\n",
    "\\begin{pmatrix}x_{t+1}\\cr \\mu_{t+1}\\cr\\end{pmatrix}\\ = M\\ \\begin{pmatrix}x_t\\cr\\mu_t\\cr\\end{pmatrix} \\tag{57.8}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "\n",
    "<a id='equation-mdefn'></a>\n",
    "$$\n",
    "M\\equiv L^{-1} N = \\begin{pmatrix}A+B Q^{-1} B^\\prime A^{\\prime-1}R &\n",
    "-B Q^{-1} B^\\prime A^{\\prime-1}\\cr -A^{\\prime -1} R & A^{\\prime -1}\\cr\\end{pmatrix}. \\tag{57.9}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afd9d16c",
   "metadata": {},
   "source": [
    "## State-Costate Dynamics\n",
    "\n",
    "We seek to solve the difference equation system  [(57.8)](#equation-eq4orig) for a sequence $ \\{x_t\\}_{t=0}^\\infty $\n",
    "that satisfies\n",
    "\n",
    "- an initial condition for $ x_0 $  \n",
    "- a terminal condition $ \\lim_{t \\rightarrow +\\infty} x_t =0 $  \n",
    "\n",
    "\n",
    "This  terminal condition reflects our desire for a **stable** solution, one that does not diverge as $ t \\rightarrow \\infty $.\n",
    "\n",
    "We inherit our wish for stability of the $ \\{x_t\\} $ sequence from a desire to maximize\n",
    "\n",
    "$$\n",
    "-\\sum_{t=0}^\\infty \\bigl[ x_t ' R x_t + u_t' Q u_t \\bigr],\n",
    "$$\n",
    "\n",
    "which requires that $ x_t' R x_t $ converge to zero as $ t \\rightarrow + \\infty $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e959fe9",
   "metadata": {},
   "source": [
    "## Reciprocal Pairs Property\n",
    "\n",
    "To proceed, we study properties of the $ (2n \\times 2n) $ matrix $ M $ defined in [(57.9)](#equation-mdefn).\n",
    "\n",
    "It helps to introduce a $ (2n \\times 2n) $ matrix\n",
    "\n",
    "$$\n",
    "J = \\begin{pmatrix}0 & -I_n\\cr I_n & 0\\cr\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "The rank of $ J $ is $ 2n $.\n",
    "\n",
    "**Definition:**  A matrix $ M $ is called **symplectic** if\n",
    "\n",
    "\n",
    "<a id='equation-lag-lqdp-eq3'></a>\n",
    "$$\n",
    "MJM^\\prime = J. \\tag{57.10}\n",
    "$$\n",
    "\n",
    "Salient properties of symplectic matrices that are readily verified include:\n",
    "\n",
    "- If $ M $ is symplectic, then $ M^2 $ is symplectic  \n",
    "- The determinant of a symplectic, then $ \\textrm{det}(M) = 1 $  \n",
    "\n",
    "\n",
    "It can be verified directly that $ M $ in equation [(57.9)](#equation-mdefn) is symplectic.\n",
    "\n",
    "It follows from equation [(57.10)](#equation-lag-lqdp-eq3) and from the fact $ J^{-1} = J^\\prime = -J $ that for any symplectic\n",
    "matrix $ M $,\n",
    "\n",
    "\n",
    "<a id='equation-lag-lqdp-eq4'></a>\n",
    "$$\n",
    "M^\\prime = J^{-1} M^{-1} J. \\tag{57.11}\n",
    "$$\n",
    "\n",
    "Equation [(57.11)](#equation-lag-lqdp-eq4) states that $ M^\\prime $ is related to the inverse of $ M $\n",
    "by a **similarity transformation**.\n",
    "\n",
    "For square matrices, recall that\n",
    "\n",
    "- similar matrices share eigenvalues  \n",
    "- eigenvalues of the inverse of a matrix are  inverses of  eigenvalues of the matrix  \n",
    "- a matrix and its transpose share eigenvalues  \n",
    "\n",
    "\n",
    "It then follows from equation [(57.11)](#equation-lag-lqdp-eq4)  that\n",
    "the eigenvalues of $ M $ occur in reciprocal pairs: if $ \\lambda $ is an\n",
    "eigenvalue of $ M $, so is $ \\lambda^{-1} $.\n",
    "\n",
    "Write equation [(57.8)](#equation-eq4orig) as\n",
    "\n",
    "\n",
    "<a id='equation-eq658'></a>\n",
    "$$\n",
    "y_{t+1} = M y_t \\tag{57.12}\n",
    "$$\n",
    "\n",
    "where $ y_t = \\begin{pmatrix}x_t\\cr \\mu_t\\cr\\end{pmatrix} $.\n",
    "\n",
    "Consider a **triangularization** of $ M $\n",
    "\n",
    "\n",
    "<a id='equation-eqn-triangledecomp'></a>\n",
    "$$\n",
    "V^{-1} M V= \\begin{pmatrix}W_{11} & W_{12} \\cr 0 & W_{22}\\cr\\end{pmatrix} \\tag{57.13}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "- each block on the right side is $ (n\\times n) $  \n",
    "- $ V $ is nonsingular  \n",
    "- all eigenvalues of $ W_{22} $ exceed $ 1 $ in modulus  \n",
    "- all  eigenvalues of $ W_{11} $ are  less than $ 1 $ in modulus  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3ae5d3f",
   "metadata": {},
   "source": [
    "## Schur decomposition\n",
    "\n",
    "The **Schur decomposition** and the **eigenvalue decomposition**\n",
    "are two  decompositions of the form [(57.13)](#equation-eqn-triangledecomp).\n",
    "\n",
    "Write equation [(57.12)](#equation-eq658)  as\n",
    "\n",
    "\n",
    "<a id='equation-eq659'></a>\n",
    "$$\n",
    "y_{t+1} = V W V^{-1} y_t. \\tag{57.14}\n",
    "$$\n",
    "\n",
    "A solution of equation [(57.14)](#equation-eq659)  for arbitrary initial condition $ y_0 $ is\n",
    "evidently\n",
    "\n",
    "\n",
    "<a id='equation-eq6510'></a>\n",
    "$$\n",
    "y_{t} = V \\left[\\begin{matrix}W^t_{11} & W_{12,t}\\cr 0 & W^t_{22}\\cr\\end{matrix}\\right]\n",
    "\\ V^{-1} y_0 \\tag{57.15}\n",
    "$$\n",
    "\n",
    "where $ W_{12,t} = W_{12} $ for $ t=1 $ and  for $ t \\geq 2 $ obeys the recursion\n",
    "\n",
    "$$\n",
    "W_{12, t} = W^{t-1}_{11} W_{12,t-1} + W_{12,t-1} W^{t-1}_{22}\n",
    "$$\n",
    "\n",
    "and where $ W^t_{ii} $ is $ W_{ii} $ raised to the $ t $th  power.\n",
    "\n",
    "Write equation [(57.15)](#equation-eq6510) as\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}y^\\ast_{1t}\\cr y^\\ast_{2t}\\cr\\end{pmatrix}\\ =\\ \\left[\\begin{matrix} W^t_{11} &\n",
    "W_{12, t}\\cr 0 & W^t_{22}\\cr\\end{matrix}\\right]\\quad \\begin{pmatrix}y^\\ast_{10}\\cr\n",
    "y^\\ast_{20}\\cr\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "where $ y^\\ast_t = V^{-1} y_t $, and in particular where\n",
    "\n",
    "\n",
    "<a id='equation-eq6511'></a>\n",
    "$$\n",
    "y^\\ast_{2t} = V^{21} x_t + V^{22} \\mu_t, \\tag{57.16}\n",
    "$$\n",
    "\n",
    "and where $ V^{ij} $ denotes the $ (i,j) $ piece of\n",
    "the partitioned $ V^{-1} $ matrix.\n",
    "\n",
    "Because $ W_{22} $ is an unstable matrix, $ y^\\ast_t $ will diverge unless $ y^\\ast_{20} = 0 $.\n",
    "\n",
    "Let $ V^{ij} $ denote the $ (i,j) $ piece of the partitioned $ V^{-1} $ matrix.\n",
    "\n",
    "To attain stability, we must impose $ y^\\ast_{20} =0 $, which from equation [(57.16)](#equation-eq6511)  implies\n",
    "\n",
    "$$\n",
    "V^{21} x_0 + V^{22} \\mu_0 = 0\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "$$\n",
    "\\mu_0 = - (V^{22})^{-1} V^{21} x_0.\n",
    "$$\n",
    "\n",
    "This equation replicates itself over\n",
    "time in the sense that it implies\n",
    "\n",
    "$$\n",
    "\\mu_t = - (V^{22})^{-1} V^{21} x_t.\n",
    "$$\n",
    "\n",
    "But notice that because $ (V^{21}\\ V^{22}) $ is the second row block of\n",
    "the inverse of $ V, $ it follows that\n",
    "\n",
    "$$\n",
    "(V^{21} \\ V^{22})\\quad \\begin{pmatrix}V_{11}\\cr V_{21}\\cr\\end{pmatrix} = 0\n",
    "$$\n",
    "\n",
    "which implies\n",
    "\n",
    "$$\n",
    "V^{21} V_{11} + V^{22} V_{21} = 0.\n",
    "$$\n",
    "\n",
    "Therefore,\n",
    "\n",
    "$$\n",
    "-(V^{22})^{-1} V^{21} = V_{21} V^{-1}_{11}.\n",
    "$$\n",
    "\n",
    "So we can write\n",
    "\n",
    "$$\n",
    "\\mu_0 = V_{21} V_{11}^{-1} x_0\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "\\mu_t = V_{21} V^{-1}_{11} x_t.\n",
    "$$\n",
    "\n",
    "However, we know  that $ \\mu_t = P x_t $,\n",
    "where $ P $ occurs in the matrix that solves the Riccati equation.\n",
    "\n",
    "Thus, the preceding argument establishes that\n",
    "\n",
    "\n",
    "<a id='equation-eqn-pvaughn'></a>\n",
    "$$\n",
    "P = V_{21} V_{11}^{-1}. \\tag{57.17}\n",
    "$$\n",
    "\n",
    "Remarkably, formula [(57.17)](#equation-eqn-pvaughn) provides us with a computationally\n",
    "efficient way of computing the positive definite  matrix $ P $ that solves the algebraic Riccati equation [(57.2)](#equation-riccati) that emerges\n",
    "from dynamic programming.\n",
    "\n",
    "This same method can be applied to compute the solution of\n",
    "any system of the form [(57.8)](#equation-eq4orig) if a solution exists, even\n",
    "if eigenvalues of $ M $ fail to occur in reciprocal pairs.\n",
    "\n",
    "The method\n",
    "will typically work so long as the eigenvalues of $ M $ split   half\n",
    "inside and half outside the unit circle.\n",
    "\n",
    "Systems in which  eigenvalues (properly adjusted for discounting) fail\n",
    "to occur in reciprocal pairs arise when the system being solved\n",
    "is an equilibrium of a model in which there are distortions that\n",
    "prevent there being any optimum problem that the equilibrium\n",
    "solves. See [[Ljungqvist and Sargent, 2018](https://python.quantecon.org/zreferences.html#id186)],  ch 12."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e61275b1",
   "metadata": {},
   "source": [
    "## Application\n",
    "\n",
    "Here we demonstrate the computation with an example which is the deterministic version of an example borrowed from this [quantecon lecture](https://python.quantecon.org/lqcontrol.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9775a54",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Model parameters\n",
    "r = 0.05\n",
    "c_bar = 2\n",
    "μ = 1\n",
    "\n",
    "# Formulate as an LQ problem\n",
    "Q = np.array([[1]])\n",
    "R = np.zeros((2, 2))\n",
    "A = [[1 + r, -c_bar + μ],\n",
    "     [0,              1]]\n",
    "B = [[-1],\n",
    "     [0]]\n",
    "\n",
    "# Construct an LQ instance\n",
    "lq = LQ(Q, R, A, B)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29b9ddf4",
   "metadata": {},
   "source": [
    "Given matrices $ A $, $ B $, $ Q $, $ R $, we can then compute $ L $, $ N $, and $ M=L^{-1}N $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d38082e4",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def construct_LNM(A, B, Q, R):\n",
    "\n",
    "    n, k = lq.n, lq.k\n",
    "\n",
    "    # construct L and N\n",
    "    L = np.zeros((2*n, 2*n))\n",
    "    L[:n, :n] = np.eye(n)\n",
    "    L[:n, n:] = B @ np.linalg.inv(Q) @ B.T\n",
    "    L[n:, n:] = A.T\n",
    "\n",
    "    N = np.zeros((2*n, 2*n))\n",
    "    N[:n, :n] = A\n",
    "    N[n:, :n] = -R\n",
    "    N[n:, n:] = np.eye(n)\n",
    "\n",
    "    # compute M\n",
    "    M = np.linalg.inv(L) @ N\n",
    "\n",
    "    return L, N, M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "760e6d11",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "L, N, M = construct_LNM(lq.A, lq.B, lq.Q, lq.R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2a3ba26",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "M"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2503e6f",
   "metadata": {},
   "source": [
    "Let’s verify that $ M $ is symplectic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb84678d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "n = lq.n\n",
    "J = np.zeros((2*n, 2*n))\n",
    "J[n:, :n] = np.eye(n)\n",
    "J[:n, n:] = -np.eye(n)\n",
    "\n",
    "M @ J @ M.T - J"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d131d8bc",
   "metadata": {},
   "source": [
    "We can compute the eigenvalues of $ M $ using `np.linalg.eigvals`, arranged in ascending order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0180d2b4",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "eigvals = sorted(np.linalg.eigvals(M))\n",
    "eigvals"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cfb9a73",
   "metadata": {},
   "source": [
    "When we apply Schur decomposition such that $ M=V W V^{-1} $, we want\n",
    "\n",
    "- the upper left block of $ W $, $ W_{11} $, to have all of its eigenvalues   less than 1 in modulus, and  \n",
    "- the lower right block $ W_{22} $ to have  eigenvalues that  exceed 1 in modulus.  \n",
    "\n",
    "\n",
    "To get what we want, let’s define a sorting function that tells `scipy.schur` to sort the corresponding eigenvalues with modulus smaller than 1 to the upper left."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f808298",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "stable_eigvals = eigvals[:n]\n",
    "\n",
    "def sort_fun(x):\n",
    "    \"Sort the eigenvalues with modules smaller than 1 to the top-left.\"\n",
    "\n",
    "    if x in stable_eigvals:\n",
    "        stable_eigvals.pop(stable_eigvals.index(x))\n",
    "        return True\n",
    "    else:\n",
    "        return False\n",
    "\n",
    "W, V, _ = schur(M, sort=sort_fun)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c6173e5",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "W"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0c7155c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "V"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86701754",
   "metadata": {},
   "source": [
    "We can check the modulus of eigenvalues of $ W_{11} $ and $ W_{22} $.\n",
    "\n",
    "Since they are both triangular matrices,  eigenvalues are the  diagonal elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2f3b73b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# W11\n",
    "np.diag(W[:n, :n])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b648ee6",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# W22\n",
    "np.diag(W[n:, n:])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8e7549a",
   "metadata": {},
   "source": [
    "The following functions wrap  $ M $ matrix construction, Schur decomposition, and stability-imposing computation of $ P $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52fbc5ee",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def stable_solution(M, verbose=True):\n",
    "    \"\"\"\n",
    "    Given a system of linear difference equations\n",
    "\n",
    "        y' = |a b| y\n",
    "        x' = |c d| x\n",
    "\n",
    "    which is potentially unstable, find the solution\n",
    "    by imposing stability.\n",
    "\n",
    "    Parameter\n",
    "    ---------\n",
    "    M : np.ndarray(float)\n",
    "        The matrix represents the linear difference equations system.\n",
    "    \"\"\"\n",
    "    n = M.shape[0] // 2\n",
    "    stable_eigvals = list(sorted(np.linalg.eigvals(M))[:n])\n",
    "\n",
    "    def sort_fun(x):\n",
    "        \"Sort the eigenvalues with modules smaller than 1 to the top-left.\"\n",
    "\n",
    "        if x in stable_eigvals:\n",
    "            stable_eigvals.pop(stable_eigvals.index(x))\n",
    "            return True\n",
    "        else:\n",
    "            return False\n",
    "\n",
    "    W, V, _ = schur(M, sort=sort_fun)\n",
    "    if verbose:\n",
    "        print('eigenvalues:\\n')\n",
    "        print('    W11: {}'.format(np.diag(W[:n, :n])))\n",
    "        print('    W22: {}'.format(np.diag(W[n:, n:])))\n",
    "\n",
    "    # compute V21 V11^{-1}\n",
    "    P = V[n:, :n] @ np.linalg.inv(V[:n, :n])\n",
    "\n",
    "    return W, V, P\n",
    "\n",
    "def stationary_P(lq, verbose=True):\n",
    "    \"\"\"\n",
    "    Computes the matrix :math:`P` that represent the value function\n",
    "\n",
    "         V(x) = x' P x\n",
    "\n",
    "    in the infinite horizon case. Computation is via imposing stability\n",
    "    on the solution path and using Schur decomposition.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    lq : qe.LQ\n",
    "        QuantEcon class for analyzing linear quadratic optimal control\n",
    "        problems of infinite horizon form.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    P : array_like(float)\n",
    "        P matrix in the value function representation.\n",
    "    \"\"\"\n",
    "\n",
    "    Q = lq.Q\n",
    "    R = lq.R\n",
    "    A = lq.A * lq.beta ** (1/2)\n",
    "    B = lq.B * lq.beta ** (1/2)\n",
    "\n",
    "    n, k = lq.n, lq.k\n",
    "\n",
    "    L, N, M = construct_LNM(A, B, Q, R)\n",
    "    W, V, P = stable_solution(M, verbose=verbose)\n",
    "\n",
    "    return P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34a40d62",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# compute P\n",
    "stationary_P(lq)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aea85e69",
   "metadata": {},
   "source": [
    "Note that the matrix $ P $ computed in this way is close to what we get from the routine in quantecon that solves an algebraic  Riccati equation by iterating to convergence on a Riccati difference equation.\n",
    "\n",
    "The small difference comes from computational errors and will decrease as we increase the maximum number of iterations or decrease the tolerance for convergence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "070d050f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "lq.stationary_values()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ceed1fb0",
   "metadata": {},
   "source": [
    "Using a Schur decomposition is much more efficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d224cd5",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%%timeit\n",
    "stationary_P(lq, verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e75dee5",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%%timeit\n",
    "lq.stationary_values()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3a5dac4",
   "metadata": {},
   "source": [
    "## Other Applications\n",
    "\n",
    "The preceding approach to imposing stability on a system  of potentially unstable linear difference equations is not limited to  linear quadratic dynamic optimization problems.\n",
    "\n",
    "For example, the same method is used in our [Stability in Linear Rational Expectations Models](https://python.quantecon.org/re_with_feedback.html#another-perspective) lecture.\n",
    "\n",
    "Let’s try to solve the model described in that lecture by applying the `stable_solution` function defined in this lecture above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92e6435d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def construct_H(ρ, λ, δ):\n",
    "    \"contruct matrix H given parameters.\"\n",
    "\n",
    "    H = np.empty((2, 2))\n",
    "    H[0, :] = ρ,δ\n",
    "    H[1, :] = - (1 - λ) / λ, 1 / λ\n",
    "\n",
    "    return H\n",
    "\n",
    "H = construct_H(ρ=.9, λ=.5, δ=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b078749b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "W, V, P = stable_solution(H)\n",
    "P"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "febeafd7",
   "metadata": {},
   "source": [
    "## Discounted Problems"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b305621",
   "metadata": {},
   "source": [
    "### Transforming States and Controls to Eliminate Discounting\n",
    "\n",
    "A pair of useful transformations allows us to convert a discounted problem into an undiscounted one.\n",
    "\n",
    "Thus, suppose that we have a discounted problem with objective\n",
    "\n",
    "$$\n",
    "- \\sum^\\infty_{t=0} \\beta^t \\biggl\\{ x^\\prime_t R x_t + u_t^\\prime Q u_t \\biggr\\}\n",
    "$$\n",
    "\n",
    "and that the state transition equation\n",
    "is again $ x_{t +1 }=Ax_t+Bu_t $.\n",
    "\n",
    "Define the transformed state and control variables\n",
    "\n",
    "- $ \\hat x_t = \\beta^{\\frac{t}{2}} x_t $  \n",
    "- $ \\hat u_t = \\beta^{\\frac{t}{2}} u_t $  \n",
    "\n",
    "\n",
    "and the transformed transition equation\n",
    "matrices\n",
    "\n",
    "- $ \\hat A = \\beta^{\\frac{1}{2}} A $  \n",
    "- $ \\hat B =  \\beta^{\\frac{1}{2}} B $  \n",
    "\n",
    "\n",
    "so that the adjusted state and control variables\n",
    "obey the transition law\n",
    "\n",
    "$$\n",
    "\\hat x_{t+1} = \\hat A \\hat x_t + \\hat B \\hat u_t.\n",
    "$$\n",
    "\n",
    "Then a discounted optimal control problem\n",
    "defined by $ A, B, R, Q, \\beta $ having  optimal policy characterized by $ P, F $ is associated with an equivalent\n",
    "undiscounted problem defined by $ \\hat A, \\hat B, Q, R $ having  optimal policy characterized by $ \\hat F, \\hat P $ that satisfy\n",
    "the following   equations:\n",
    "\n",
    "$$\n",
    "\\hat F=(Q+B'\\hat PB)^{-1}\\hat B'P \\hat A\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "\\hat P=R+\\hat A'P \\hat A-\\hat A'P \\hat B(Q+B'\\hat P \\hat B)^{-1} \\hat B'P \\hat A\n",
    "$$\n",
    "\n",
    "It follows immediately from the definitions of $ \\hat A, \\hat B $ that $ \\hat F = F $ and $ \\hat P = P $.\n",
    "\n",
    "By exploiting these transformations,  we can solve a discounted problem by solving an associated undiscounted problem.\n",
    "\n",
    "In particular, we can first transform a discounted LQ problem to an undiscounted one and then solve that  discounted optimal regulator problem using the Lagrangian and invariant subspace methods described above.\n",
    "\n",
    "For example, when $ \\beta=\\frac{1}{1+r} $, we can solve for $ P $ with $ \\hat{A}=\\beta^{1/2} A $ and $ \\hat{B}=\\beta^{1/2} B $.\n",
    "\n",
    "These settings are adopted by default in the function `stationary_P` defined above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "064dc7d7",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "β = 1 / (1 + r)\n",
    "lq.beta = β"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39becdd4",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "stationary_P(lq)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf1a2f9f",
   "metadata": {},
   "source": [
    "We can verify that the solution agrees with one that comes from applying the  routine `LQ.stationary_values` in the  quantecon package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0f02337",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "lq.stationary_values()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "796a2c24",
   "metadata": {},
   "source": [
    "### Lagrangian for Discounted Problem\n",
    "\n",
    "For several purposes, it is useful  explicitly briefly to describe\n",
    "a Lagrangian for a discounted problem.\n",
    "\n",
    "Thus, for the discounted optimal linear regulator problem,\n",
    "form the Lagrangian\n",
    "\n",
    "\n",
    "<a id='equation-eq661'></a>\n",
    "$$\n",
    "{\\cal{L}} = - \\sum^\\infty_{t=0} \\beta^t \\biggl\\{ x^\\prime_t R x_t + u_t^\\prime Q u_t\n",
    "+ 2 \\beta \\mu^\\prime_{t+1} [A x_t + B u_t - x_{t+1}]\\biggr\\} \\tag{57.18}\n",
    "$$\n",
    "\n",
    "where $ 2 \\mu_{t+1} $ is a vector of Lagrange multipliers on the state vector $ x_{t+1} $.\n",
    "\n",
    "First-order conditions for maximization with respect\n",
    "to $ \\{u_t,x_{t+1}\\}_{t=0}^\\infty $ are\n",
    "\n",
    "\n",
    "<a id='equation-eq662'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "2 Q u_t &+ 2  \\beta B^\\prime \\mu_{t+1} = 0 \\ ,\\ t \\geq 0 \\cr \\mu_t &= R x_t + \\beta A^\\prime \\mu_{t+1}\\ ,\\ t\\geq 1.\\cr\n",
    "\\end{aligned} \\tag{57.19}\n",
    "$$\n",
    "\n",
    "Define $ 2 \\mu_0 $ to be the vector of shadow prices of $ x_0 $ and apply an envelope condition to\n",
    "[(57.18)](#equation-eq661) to  deduce that\n",
    "\n",
    "$$\n",
    "\\mu_0 = R x_0 + \\beta A' \\mu_1 ,\n",
    "$$\n",
    "\n",
    "which is a time $ t=0 $ counterpart to the second equation of system [(57.19)](#equation-eq662).\n",
    "\n",
    "Proceeding as we did above with  the undiscounted system  [(57.5)](#equation-lag-lqdp-eq2), we can rearrange the first-order conditions into the\n",
    "system\n",
    "\n",
    "\n",
    "<a id='equation-eq663'></a>\n",
    "$$\n",
    "\\left[\\begin{matrix} I & \\beta B Q^{-1} B' \\cr\n",
    "             0 & \\beta A' \\end{matrix}\\right]\n",
    "\\left[\\begin{matrix} x_{t+1} \\cr \\mu_{t+1} \\end{matrix}\\right] =\n",
    "\\left[\\begin{matrix} A & 0 \\cr\n",
    "             - R & I \\end{matrix}\\right] \n",
    "\\left[\\begin{matrix} x_t \\cr \\mu_t \\end{matrix}\\right] \\tag{57.20}\n",
    "$$\n",
    "\n",
    "which in the special case that $ \\beta = 1 $ agrees with equation [(57.5)](#equation-lag-lqdp-eq2), as expected.\n",
    "\n",
    "By staring at system [(57.20)](#equation-eq663), we can infer  identities that shed light on the structure of optimal linear regulator problems, some of which will be useful in [this lecture](https://python-advanced.quantecon.org/dyn_stack.html) when we apply and  extend the methods of this lecture to study Stackelberg and Ramsey problems.\n",
    "\n",
    "First, note that the first block of equation system [(57.20)](#equation-eq663) asserts that when  $ \\mu_{t+1} = P x_{t+1} $, then\n",
    "\n",
    "$$\n",
    "(I + \\beta Q^{-1} B' P B P ) x_{t+1} = A x_t,\n",
    "$$\n",
    "\n",
    "which  can be rearranged to sbe\n",
    "\n",
    "$$\n",
    "x_{t+1} = (I + \\beta B Q^{-1} B' P)^{-1}  A x_t .\n",
    "$$\n",
    "\n",
    "This expression for the optimal closed loop dynamics of the state  must agree with an alternative expression that we had derived with dynamic programming, namely,\n",
    "\n",
    "$$\n",
    "x_{t+1} = (A - BF) x_t .\n",
    "$$\n",
    "\n",
    "But using\n",
    "\n",
    "\n",
    "<a id='equation-eqn-optimalfformula'></a>\n",
    "$$\n",
    "F=\\beta (Q+\\beta B'PB)^{-1} B'PA \\tag{57.21}\n",
    "$$\n",
    "\n",
    "it follows that\n",
    "\n",
    "$$\n",
    "A- B F = (I - \\beta B (Q+ \\beta B' P B)^{-1} B' P) A .\n",
    "$$\n",
    "\n",
    "Thus, our two expressions for the\n",
    "closed loop dynamics  agree if and only if\n",
    "\n",
    "\n",
    "<a id='equation-eqn-twofeedbackloops'></a>\n",
    "$$\n",
    "(I + \\beta B Q^{-1} B' P )^{-1} =    (I - \\beta B (Q+\\beta  B' P B)^{-1} B' P) . \\tag{57.22}\n",
    "$$\n",
    "\n",
    "Matrix  equation [(57.22)](#equation-eqn-twofeedbackloops) can be verified by applying a partitioned inverse formula.\n",
    "\n",
    ">**Note**\n",
    ">\n",
    ">Just use the formula $ (a - b d^{-1} c)^{-1} = a^{-1} + a^{-1} b (d - c a^{-1} b)^{-1} c a^{-1} $ for appropriate choices of the matrices $ a, b, c, d $.\n",
    "\n",
    "Next, note that for *any*  fixed $ F $ for which eigenvalues of $ A- BF $ are less than $ \\frac{1}{\\beta} $ in modulus, the value function associated with using this rule forever is $ - x_0 \\tilde P x_0 $ where  $ \\tilde P $ obeys the following  matrix  equation:\n",
    "\n",
    "\n",
    "<a id='equation-eq666'></a>\n",
    "$$\n",
    "\\tilde P = (R + F' Q F) + \\beta (A - B F)' P (A - BF) . \\tag{57.23}\n",
    "$$\n",
    "\n",
    "Evidently, $ \\tilde P = P $ only when $ F $ obeys formula [(57.21)](#equation-eqn-optimalfformula).\n",
    "\n",
    "Next, note that  the second equation of system [(57.20)](#equation-eq663) implies the “forward looking” equation for the Lagrange multiplier\n",
    "\n",
    "$$\n",
    "\\mu_t = R x_t + \\beta A' \\mu_{t+1}\n",
    "$$\n",
    "\n",
    "whose solution is\n",
    "\n",
    "$$\n",
    "\\mu_t = P x_t ,\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "\n",
    "<a id='equation-eq667'></a>\n",
    "$$\n",
    "P = R + \\beta A' P (A - BF) \\tag{57.24}\n",
    "$$\n",
    "\n",
    "where we must require that $ F $ obeys equation [(57.21)](#equation-eqn-optimalfformula).\n",
    "\n",
    "Equations [(57.23)](#equation-eq666) and [(57.24)](#equation-eq667) provide different perspectives on the optimal value function."
   ]
  }
 ],
 "metadata": {
  "date": 1731921899.9411957,
  "filename": "lagrangian_lqdp.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Lagrangian for LQ Control"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}