markov_perf.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='markov-perf'></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": [
    "# Markov Perfect Equilibrium\n",
    "\n",
    "\n",
    "<a id='index-0'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Markov Perfect Equilibrium](#Markov-Perfect-Equilibrium)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [Background](#Background)  \n",
    "  - [Linear Markov Perfect Equilibria](#Linear-Markov-Perfect-Equilibria)  \n",
    "  - [Application](#Application)  \n",
    "  - [Exercises](#Exercises)  \n",
    "  - [Solutions](#Solutions)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to what&rsquo;s in Anaconda, this lecture will need the following libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": true
   },
   "outputs": [],
   "source": [
    "!pip install --upgrade quantecon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "\n",
    "<a id='index-1'></a>\n",
    "This lecture describes the concept of Markov perfect equilibrium.\n",
    "\n",
    "Markov perfect equilibrium is a key notion for analyzing economic problems involving dynamic strategic interaction, and a cornerstone of applied game theory.\n",
    "\n",
    "In this lecture, we teach Markov perfect equilibrium by example.\n",
    "\n",
    "We will focus on settings with\n",
    "\n",
    "- two players  \n",
    "- quadratic payoff functions  \n",
    "- linear transition rules for the state  \n",
    "\n",
    "\n",
    "Other references include chapter 7 of [[LS18]](https://lectures.quantecon.org/zreferences.html#ljungqvist2012).\n",
    "\n",
    "Let&rsquo;s start with some standard imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import quantecon as qe\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Background\n",
    "\n",
    "\n",
    "<a id='index-2'></a>\n",
    "Markov perfect equilibrium is a refinement of the concept of Nash equilibrium.\n",
    "\n",
    "It is used to study settings where multiple decision-makers interact non-cooperatively over time, each seeking to pursue its own objective.\n",
    "\n",
    "The agents in the model face a common state vector, the time path of which is influenced by &ndash; and influences &ndash; their decisions.\n",
    "\n",
    "In particular, the transition law for the state that confronts each agent is affected by decision rules of other agents.\n",
    "\n",
    "Individual payoff maximization requires that each agent solve a dynamic programming problem that includes  this transition law.\n",
    "\n",
    "Markov perfect equilibrium prevails when no agent wishes to revise its policy, taking as given the policies of all other agents.\n",
    "\n",
    "Well known examples include\n",
    "\n",
    "- Choice of price, output, location or capacity for firms in an industry (e.g., [[EP95]](https://lectures.quantecon.org/zreferences.html#ericson1995markov), [[Rya12]](https://lectures.quantecon.org/zreferences.html#ryan2012costs), [[DS10]](https://lectures.quantecon.org/zreferences.html#doraszelski2010computable)).  \n",
    "- Rate of extraction from a shared natural resource, such as a fishery (e.g., [[LM80]](https://lectures.quantecon.org/zreferences.html#levhari1980great), [[VL11]](https://lectures.quantecon.org/zreferences.html#van2011dynamic)).  \n",
    "\n",
    "\n",
    "Let&rsquo;s examine a model of the first type."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: A Duopoly Model\n",
    "\n",
    "Two firms are the only producers of a good the demand for which is governed by a linear inverse demand function\n",
    "\n",
    "\n",
    "<a id='equation-game2'></a>\n",
    "$$\n",
    "p = a_0 - a_1 (q_1 +  q_2) \\tag{1}\n",
    "$$\n",
    "\n",
    "Here $ p = p_t $ is the price of the good, $ q_i = q_{it} $ is the output of firm $ i=1,2 $ at time $ t $ and $ a_0 > 0, a_1 >0 $.\n",
    "\n",
    "In [(1)](#equation-game2) and what follows,\n",
    "\n",
    "- the time subscript is suppressed when possible to simplify notation  \n",
    "- $ \\hat x $ denotes a next period value of variable $ x $  \n",
    "\n",
    "\n",
    "Each firm recognizes that its output affects total output and therefore the market price.\n",
    "\n",
    "The one-period payoff function of firm $ i $ is price times quantity minus adjustment costs:\n",
    "\n",
    "\n",
    "<a id='equation-game1'></a>\n",
    "$$\n",
    "\\pi_i = p q_i - \\gamma (\\hat q_i - q_i)^2, \\quad \\gamma > 0 , \\tag{2}\n",
    "$$\n",
    "\n",
    "Substituting the inverse demand curve [(1)](#equation-game2) into [(2)](#equation-game1) lets us express the one-period payoff as\n",
    "\n",
    "\n",
    "<a id='equation-game3'></a>\n",
    "$$\n",
    "\\pi_i(q_i, q_{-i}, \\hat q_i) = a_0 q_i - a_1 q_i^2 - a_1 q_i q_{-i} - \\gamma (\\hat q_i - q_i)^2 , \\tag{3}\n",
    "$$\n",
    "\n",
    "where $ q_{-i} $ denotes the output of the firm other than $ i $.\n",
    "\n",
    "The objective of the firm is to maximize $ \\sum_{t=0}^\\infty \\beta^t \\pi_{it} $.\n",
    "\n",
    "Firm $ i $ chooses a decision rule that sets next period quantity $ \\hat q_i $ as a function $ f_i $ of the current state $ (q_i, q_{-i}) $.\n",
    "\n",
    "An essential aspect of a Markov perfect equilibrium is that each firm takes the decision rule of the other firm as known and given.\n",
    "\n",
    "Given $ f_{-i} $, the Bellman equation of firm $ i $ is\n",
    "\n",
    "\n",
    "<a id='equation-game4'></a>\n",
    "$$\n",
    "v_i(q_i, q_{-i}) = \\max_{\\hat q_i}\n",
    "   \\left\\{\\pi_i (q_i, q_{-i}, \\hat q_i) + \\beta v_i(\\hat q_i, f_{-i}(q_{-i}, q_i)) \\right\\} \\tag{4}\n",
    "$$\n",
    "\n",
    "**Definition**  A *Markov perfect equilibrium* of the duopoly model is a pair of value functions $ (v_1, v_2) $ and a pair of policy functions $ (f_1, f_2) $ such that, for each $ i \\in \\{1, 2\\} $ and each possible state,\n",
    "\n",
    "- The value function $ v_i $ satisfies the Bellman equation [(4)](#equation-game4).  \n",
    "- The maximizer on the right side of [(4)](#equation-game4) is equal to $ f_i(q_i, q_{-i}) $.  \n",
    "\n",
    "\n",
    "The adjective &ldquo;Markov&rdquo; denotes that the equilibrium decision rules depend only on the current values of the state variables, not other parts of their histories.\n",
    "\n",
    "&ldquo;Perfect&rdquo; means complete, in the sense that the equilibrium is constructed by backward induction and hence builds in optimizing behavior for each firm at all possible future states.\n",
    "\n",
    "> - These include many states that will not be reached when we iterate forward on the pair of equilibrium strategies $ f_i $ starting from a given initial state.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computation\n",
    "\n",
    "One strategy for computing a Markov perfect equilibrium is iterating to convergence on pairs of Bellman equations and decision rules.\n",
    "\n",
    "In particular, let $ v_i^j,f_i^j $ be the value function and policy function for firm $ i $ at the $ j $-th iteration.\n",
    "\n",
    "Imagine constructing the iterates\n",
    "\n",
    "\n",
    "<a id='equation-game6'></a>\n",
    "$$\n",
    "v_i^{j+1}(q_i, q_{-i}) = \\max_{\\hat q_i}\n",
    "   \\left\\{\\pi_i (q_i, q_{-i}, \\hat q_i) + \\beta v^j_i(\\hat q_i, f_{-i}(q_{-i}, q_i)) \\right\\} \\tag{5}\n",
    "$$\n",
    "\n",
    "These iterations can be challenging to implement computationally.\n",
    "\n",
    "However, they simplify for the case in which the one-period payoff functions are quadratic and the transition laws are linear &mdash; which takes us to our next topic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Markov Perfect Equilibria\n",
    "\n",
    "\n",
    "<a id='index-3'></a>\n",
    "As we saw in the duopoly example, the study of Markov perfect equilibria in games with two players leads us to an interrelated pair of Bellman equations.\n",
    "\n",
    "In linear-quadratic dynamic games, these &ldquo;stacked Bellman equations&rdquo; become &ldquo;stacked Riccati equations&rdquo; with a tractable mathematical structure.\n",
    "\n",
    "We&rsquo;ll lay out that structure in a general setup and then apply it to some simple problems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coupled Linear Regulator Problems\n",
    "\n",
    "We consider a general linear-quadratic regulator game with two players.\n",
    "\n",
    "For convenience, we&rsquo;ll start with a finite horizon formulation, where $ t_0 $ is the initial date and $ t_1 $ is the common terminal date.\n",
    "\n",
    "Player $ i $ takes $ \\{u_{-it}\\} $ as given and minimizes\n",
    "\n",
    "\n",
    "<a id='equation-orig-1'></a>\n",
    "$$\n",
    "\\sum_{t=t_0}^{t_1 - 1}\n",
    "\\beta^{t - t_0}\n",
    "\\left\\{\n",
    "    x_t' R_i x_t +\n",
    "    u_{it}' Q_i u_{it} +\n",
    "    u_{-it}' S_i u_{-it} +\n",
    "    2 x_t' W_i u_{it} +\n",
    "    2 u_{-it}' M_i u_{it}\n",
    "\\right\\} \\tag{6}\n",
    "$$\n",
    "\n",
    "while the state evolves according to\n",
    "\n",
    "\n",
    "<a id='equation-orig-0'></a>\n",
    "$$\n",
    "x_{t+1} = A x_t + B_1 u_{1t} + B_2 u_{2t} \\tag{7}\n",
    "$$\n",
    "\n",
    "Here\n",
    "\n",
    "- $ x_t $ is an $ n \\times 1 $ state vector and  $ u_{it} $ is a $ k_i \\times 1 $ vector of controls for player $ i $  \n",
    "- $ R_i $ is $ n \\times n $  \n",
    "- $ S_i $ is $ k_{-i} \\times k_{-i} $  \n",
    "- $ Q_i $ is $ k_i \\times k_i $  \n",
    "- $ W_i $ is $ n \\times k_i $  \n",
    "- $ M_i $ is $ k_{-i} \\times k_i $  \n",
    "- $ A $ is $ n \\times n $  \n",
    "- $ B_i $ is $ n \\times k_i $  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing Equilibrium\n",
    "\n",
    "We formulate a linear Markov perfect equilibrium as follows.\n",
    "\n",
    "Player $ i $ employs linear decision rules $ u_{it} = - F_{it} x_t $, where $ F_{it} $ is a $ k_i \\times n $ matrix.\n",
    "\n",
    "A Markov perfect equilibrium is a pair of sequences $ \\{F_{1t}, F_{2t}\\} $ over $ t = t_0, \\ldots, t_1 - 1 $ such that\n",
    "\n",
    "- $ \\{F_{1t}\\} $ solves player 1&rsquo;s problem, taking $ \\{F_{2t}\\} $ as given, and  \n",
    "- $ \\{F_{2t}\\} $ solves player 2&rsquo;s problem, taking $ \\{F_{1t}\\} $ as given  \n",
    "\n",
    "\n",
    "If we take $ u_{2t} = - F_{2t} x_t $ and substitute it into [(6)](#equation-orig-1) and [(7)](#equation-orig-0), then player 1&rsquo;s problem becomes minimization of\n",
    "\n",
    "\n",
    "<a id='equation-eq-mpe-p1p'></a>\n",
    "$$\n",
    "\\sum_{t=t_0}^{t_1 - 1}\n",
    "\\beta^{t - t_0}\n",
    "    \\left\\{\n",
    "    x_t' \\Pi_{1t} x_t +\n",
    "    u_{1t}' Q_1 u_{1t} +\n",
    "    2 u_{1t}' \\Gamma_{1t} x_t\n",
    "    \\right\\} \\tag{8}\n",
    "$$\n",
    "\n",
    "subject to\n",
    "\n",
    "\n",
    "<a id='equation-eq-mpe-p1d'></a>\n",
    "$$\n",
    "x_{t+1} = \\Lambda_{1t} x_t + B_1 u_{1t}, \\tag{9}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "- $ \\Lambda_{it} := A - B_{-i} F_{-it} $  \n",
    "- $ \\Pi_{it} := R_i + F_{-it}' S_i F_{-it} $  \n",
    "- $ \\Gamma_{it} := W_i' - M_i' F_{-it} $  \n",
    "\n",
    "\n",
    "This is an LQ dynamic programming problem that can be solved by working backwards.\n",
    "\n",
    "The policy rule that solves this problem is\n",
    "\n",
    "\n",
    "<a id='equation-orig-3'></a>\n",
    "$$\n",
    "F_{1t}\n",
    "= (Q_1 + \\beta B_1' P_{1t+1} B_1)^{-1}\n",
    "(\\beta B_1' P_{1t+1} \\Lambda_{1t} + \\Gamma_{1t}) \\tag{10}\n",
    "$$\n",
    "\n",
    "where $ P_{1t} $ solves the matrix Riccati difference equation\n",
    "\n",
    "\n",
    "<a id='equation-orig-4'></a>\n",
    "$$\n",
    "P_{1t} =\n",
    "\\Pi_{1t} - (\\beta B_1' P_{1t+1} \\Lambda_{1t} +\n",
    "\\Gamma_{1t})' (Q_1 + \\beta B_1' P_{1t+1} B_1)^{-1}\n",
    "(\\beta B_1' P_{1t+1} \\Lambda_{1t} + \\Gamma_{1t}) +\n",
    "\\beta \\Lambda_{1t}' P_{1t+1} \\Lambda_{1t} \\tag{11}\n",
    "$$\n",
    "\n",
    "Similarly, the policy that solves player 2&rsquo;s problem is\n",
    "\n",
    "\n",
    "<a id='equation-orig-5'></a>\n",
    "$$\n",
    "F_{2t} = (Q_2 + \\beta B_2' P_{2t+1} B_2)^{-1}\n",
    "(\\beta B_2' P_{2t+1} \\Lambda_{2t} + \\Gamma_{2t}) \\tag{12}\n",
    "$$\n",
    "\n",
    "where $ P_{2t} $ solves\n",
    "\n",
    "\n",
    "<a id='equation-orig-6'></a>\n",
    "$$\n",
    "P_{2t} =\n",
    "\\Pi_{2t} - (\\beta B_2' P_{2t+1} \\Lambda_{2t} +\n",
    "\\Gamma_{2t})' (Q_2 + \\beta B_2' P_{2t+1} B_2)^{-1}\n",
    "(\\beta B_2' P_{2t+1} \\Lambda_{2t} + \\Gamma_{2t}) +\n",
    "\\beta \\Lambda_{2t}' P_{2t+1} \\Lambda_{2t} \\tag{13}\n",
    "$$\n",
    "\n",
    "Here in all cases $ t = t_0, \\ldots, t_1 - 1 $ and the terminal conditions are $ P_{it_1} = 0 $.\n",
    "\n",
    "The solution procedure is to use equations [(10)](#equation-orig-3), [(11)](#equation-orig-4), [(12)](#equation-orig-5), and [(13)](#equation-orig-6), and &ldquo;work backwards&rdquo; from time $ t_1 - 1 $.\n",
    "\n",
    "Since we&rsquo;re working backward, $ P_{1t+1} $ and $ P_{2t+1} $ are taken as given at each stage.\n",
    "\n",
    "Moreover, since\n",
    "\n",
    "- some terms on the right-hand side of [(10)](#equation-orig-3) contain $ F_{2t} $  \n",
    "- some terms on the right-hand side of [(12)](#equation-orig-5) contain $ F_{1t} $  \n",
    "\n",
    "\n",
    "we need to solve these $ k_1 + k_2 $ equations simultaneously."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Key Insight\n",
    "\n",
    "A key insight is that  equations  [(10)](#equation-orig-3) and [(12)](#equation-orig-5) are linear in $ F_{1t} $ and $ F_{2t} $.\n",
    "\n",
    "After these equations have been solved, we can take  $ F_{it} $ and solve for $ P_{it} $ in [(11)](#equation-orig-4) and [(13)](#equation-orig-6)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Infinite Horizon\n",
    "\n",
    "We often want to compute the solutions of such games for infinite horizons, in the hope that the decision rules $ F_{it} $ settle down to be time-invariant as $ t_1 \\rightarrow +\\infty $.\n",
    "\n",
    "In practice, we usually fix $ t_1 $ and compute the equilibrium of an infinite horizon game by driving $ t_0 \\rightarrow - \\infty $.\n",
    "\n",
    "This is the approach we adopt in the next section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation\n",
    "\n",
    "We use the function [nnash](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lqnash.py) from [QuantEcon.py](http://quantecon.org/python_index.html) that computes a Markov perfect equilibrium of the infinite horizon linear-quadratic dynamic game in the manner described above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Application\n",
    "\n",
    "\n",
    "<a id='index-4'></a>\n",
    "Let&rsquo;s use these procedures to treat some applications, starting with the duopoly model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A Duopoly Model\n",
    "\n",
    "To map the duopoly model into  coupled linear-quadratic dynamic programming problems, define the state\n",
    "and controls as\n",
    "\n",
    "$$\n",
    "x_t :=\n",
    "\\begin{bmatrix}\n",
    "    1 \\\\\n",
    "    q_{1t} \\\\\n",
    "    q_{2t}\n",
    "\\end{bmatrix}\n",
    "\\quad \\text{and} \\quad\n",
    "u_{it} :=\n",
    "q_{i,t+1} - q_{it}, \\quad i=1,2\n",
    "$$\n",
    "\n",
    "If we write\n",
    "\n",
    "$$\n",
    "x_t' R_i x_t + u_{it}' Q_i u_{it}\n",
    "$$\n",
    "\n",
    "where $ Q_1 = Q_2 = \\gamma $,\n",
    "\n",
    "$$\n",
    "R_1 :=\n",
    "\\begin{bmatrix}\n",
    "   0              & -\\frac{a_0}{2}  & 0 \\\\\n",
    "   -\\frac{a_0}{2} &  a_1            &  \\frac{a_1}{2} \\\\\n",
    "   0              &   \\frac{a_1}{2} & 0\n",
    "\\end{bmatrix}\n",
    "\\quad \\text{and} \\quad\n",
    "R_2 :=\n",
    "\\begin{bmatrix}\n",
    "   0              & 0             & -\\frac{a_0}{2} \\\\\n",
    "   0              & 0             & \\frac{a_1}{2} \\\\\n",
    "   -\\frac{a_0}{2} & \\frac{a_1}{2} & a_1\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "then we recover the  one-period  payoffs in expression [(3)](#equation-game3).\n",
    "\n",
    "The law of motion for the state $ x_t $ is $ x_{t+1} = A x_t + B_1 u_{1t} + B_2 u_{2t} $ where\n",
    "\n",
    "$$\n",
    "A :=\n",
    "\\begin{bmatrix}\n",
    "   1 & 0 & 0 \\\\\n",
    "   0 & 1 & 0 \\\\\n",
    "   0 & 0 & 1\n",
    "\\end{bmatrix},\n",
    "\\quad\n",
    "B_1 :=\n",
    "\\begin{bmatrix}\n",
    "       0 \\\\\n",
    "       1 \\\\\n",
    "       0\n",
    "\\end{bmatrix},\n",
    "\\quad\n",
    "B_2 :=\n",
    "\\begin{bmatrix}\n",
    "            0 \\\\\n",
    "            0 \\\\\n",
    "            1\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "The optimal decision rule of firm $ i $ will take the form $ u_{it} = - F_i x_t $, inducing the following closed-loop system for the evolution of $ x $ in the Markov perfect equilibrium:\n",
    "\n",
    "\n",
    "<a id='equation-eq-mpe-cle'></a>\n",
    "$$\n",
    "x_{t+1} = (A - B_1 F_1 -B_1 F_2 ) x_t \\tag{14}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameters and Solution\n",
    "\n",
    "Consider the previously presented duopoly model with parameter values of:\n",
    "\n",
    "- $ a_0 = 10 $  \n",
    "- $ a_1 = 2 $  \n",
    "- $ \\beta = 0.96 $  \n",
    "- $ \\gamma = 12 $  \n",
    "\n",
    "\n",
    "From these, we compute the infinite horizon MPE using the preceding code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "@authors: Chase Coleman, Thomas Sargent, John Stachurski\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "import quantecon as qe\n",
    "\n",
    "# Parameters\n",
    "a0 = 10.0\n",
    "a1 = 2.0\n",
    "&beta; = 0.96\n",
    "&gamma; = 12.0\n",
    "\n",
    "# In LQ form\n",
    "A = np.eye(3)\n",
    "B1 = np.array([[0.], [1.], [0.]])\n",
    "B2 = np.array([[0.], [0.], [1.]])\n",
    "\n",
    "\n",
    "R1 = [[      0.,     -a0 / 2,          0.],\n",
    "      [-a0 / 2.,          a1,     a1 / 2.],\n",
    "      [       0,     a1 / 2.,          0.]]\n",
    "\n",
    "R2 = [[     0.,           0.,      -a0 / 2],\n",
    "      [     0.,           0.,      a1 / 2.],\n",
    "      [-a0 / 2,      a1 / 2.,           a1]]\n",
    "\n",
    "Q1 = Q2 = &gamma;\n",
    "S1 = S2 = W1 = W2 = M1 = M2 = 0.0\n",
    "\n",
    "# Solve using QE's nnash function\n",
    "F1, F2, P1, P2 = qe.nnash(A, B1, B2, R1, R2, Q1, \n",
    "                          Q2, S1, S2, W1, W2, M1, \n",
    "                          M2, beta=&beta;)\n",
    "\n",
    "# Display policies\n",
    "print(\"Computed policies for firm 1 and firm 2:\\n\")\n",
    "print(f\"F1 = {F1}\")\n",
    "print(f\"F2 = {F2}\")\n",
    "print(\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running the code produces the following output.\n",
    "\n",
    "One way to see that $ F_i $ is indeed optimal for firm $ i $ taking $ F_2 $ as given is to use [QuantEcon.py](http://quantecon.org/python_index.html)&rsquo;s LQ class.\n",
    "\n",
    "In particular, let&rsquo;s take F2 as computed above, plug it into [(8)](#equation-eq-mpe-p1p) and [(9)](#equation-eq-mpe-p1d) to get firm 1&rsquo;s problem and solve it using LQ.\n",
    "\n",
    "We hope that the resulting policy will agree with F1 as computed above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "&Lambda;1 = A - B2 @ F2\n",
    "lq1 = qe.LQ(Q1, R1, &Lambda;1, B1, beta=&beta;)\n",
    "P1_ih, F1_ih, d = lq1.stationary_values()\n",
    "F1_ih"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is close enough for rock and roll, as they say in the trade.\n",
    "\n",
    "Indeed, np.allclose agrees with our assessment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "np.allclose(F1, F1_ih)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dynamics\n",
    "\n",
    "Let&rsquo;s now investigate the dynamics of price and output in this simple duopoly model under the MPE policies.\n",
    "\n",
    "Given our optimal policies $ F1 $ and $ F2 $, the state evolves according to [(14)](#equation-eq-mpe-cle).\n",
    "\n",
    "The following program\n",
    "\n",
    "- imports $ F1 $ and $ F2 $ from the previous program along with all parameters.  \n",
    "- computes the evolution of $ x_t $ using [(14)](#equation-eq-mpe-cle).  \n",
    "- extracts and plots industry output $ q_t = q_{1t} + q_{2t} $ and price $ p_t = a_0 - a_1 q_t $.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "AF = A - B1 @ F1 - B2 @ F2\n",
    "n = 20\n",
    "x = np.empty((3, n))\n",
    "x[:, 0] = 1, 1, 1\n",
    "for t in range(n-1):\n",
    "    x[:, t+1] = AF @ x[:, t]\n",
    "q1 = x[1, :]\n",
    "q2 = x[2, :]\n",
    "q = q1 + q2       # Total output, MPE\n",
    "p = a0 - a1 * q   # Price, MPE\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(9, 5.8))\n",
    "ax.plot(q, 'b-', lw=2, alpha=0.75, label='total output')\n",
    "ax.plot(p, 'g-', lw=2, alpha=0.75, label='price')\n",
    "ax.set_title('Output and prices, duopoly MPE')\n",
    "ax.legend(frameon=False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the initial condition has been set to $ q_{10} = q_{20} = 1.0 $.\n",
    "\n",
    "To gain some perspective we can compare this to what happens in the monopoly case.\n",
    "\n",
    "The first panel in the next figure compares output of the monopolist and industry output under the MPE, as a function of time.\n",
    "\n",
    "The second panel shows analogous curves for price.\n",
    "\n",
    "\n",
    "<a id='mpe-vs-monopolist'></a>\n",
    "<img src=\"https://s3-ap-southeast-2.amazonaws.com/lectures.quantecon.org/py/_static/lecture_specific/markov_perf/mpe_vs_monopolist.png\" style=\"\">\n",
    "\n",
    "  \n",
    "Here parameters are the same as above for both the MPE and monopoly solutions.\n",
    "\n",
    "The monopolist initial condition is $ q_0 = 2.0 $ to mimic the industry initial condition $ q_{10} = q_{20} = 1.0 $ in the MPE case.\n",
    "\n",
    "As expected, output is higher and prices are lower under duopoly than monopoly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Replicate the [pair of figures](#mpe-vs-monopolist) showing the comparison of output and prices for the monopolist and duopoly under MPE.\n",
    "\n",
    "Parameters are as in duopoly_mpe.py and you can use that code to compute MPE policies under duopoly.\n",
    "\n",
    "The optimal policy in the monopolist case can be computed using [QuantEcon.py](http://quantecon.org/python_index.html)&rsquo;s LQ class."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "\n",
    "In this exercise, we consider a slightly more sophisticated duopoly problem.\n",
    "\n",
    "It takes the form of infinite horizon linear-quadratic game proposed by Judd [[Jud90]](https://lectures.quantecon.org/zreferences.html#judd1990).\n",
    "\n",
    "Two firms set prices and quantities of two goods interrelated through their demand curves.\n",
    "\n",
    "Relevant variables are defined as follows:\n",
    "\n",
    "- $ I_{it} $ = inventories of firm $ i $ at beginning of $ t $  \n",
    "- $ q_{it} $ = production of firm $ i $ during period $ t $  \n",
    "- $ p_{it} $ = price charged by firm $ i $ during period $ t $  \n",
    "- $ S_{it} $ = sales made by firm $ i $ during period $ t $  \n",
    "- $ E_{it} $ = costs of production of firm $ i $ during period $ t $  \n",
    "- $ C_{it} $ = costs of carrying inventories for firm $ i $ during $ t $  \n",
    "\n",
    "\n",
    "The firms&rsquo; cost functions are\n",
    "\n",
    "- $ C_{it} = c_{i1} + c_{i2} I_{it} + 0.5 c_{i3} I_{it}^2 $  \n",
    "- $ E_{it} = e_{i1} + e_{i2}q_{it} + 0.5 e_{i3} q_{it}^2 $ where $ e_{ij}, c_{ij} $ are positive scalars  \n",
    "\n",
    "\n",
    "Inventories obey the laws of motion\n",
    "\n",
    "$$\n",
    "I_{i,t+1} = (1 - \\delta)  I_{it} + q_{it} - S_{it}\n",
    "$$\n",
    "\n",
    "Demand is governed by the linear schedule\n",
    "\n",
    "$$\n",
    "S_t = D p_{it} + b\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "- $ S_t = \\begin{bmatrix} S_{1t} & S_{2t} \\end{bmatrix}' $  \n",
    "- $ D $ is a $ 2\\times 2 $ negative definite matrix and  \n",
    "- $ b $ is a vector of constants  \n",
    "\n",
    "\n",
    "Firm $ i $ maximizes the undiscounted sum\n",
    "\n",
    "$$\n",
    "\\lim_{T \\to \\infty}\\ {1 \\over T}\\   \\sum^T_{t=0}\\   \\left( p_{it} S_{it} - E_{it} - C_{it} \\right)\n",
    "$$\n",
    "\n",
    "We can convert this to a linear-quadratic problem by taking\n",
    "\n",
    "$$\n",
    "u_{it} =\n",
    "\\begin{bmatrix}\n",
    "    p_{it} \\\\\n",
    "    q_{it}\n",
    "\\end{bmatrix}\n",
    "\\quad \\text{and} \\quad\n",
    "x_t =\n",
    "\\begin{bmatrix}\n",
    "    I_{1t} \\\\\n",
    "    I_{2t} \\\\\n",
    "    1\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Decision rules for price and quantity take the form $ u_{it} = -F_i  x_t $.\n",
    "\n",
    "The Markov perfect equilibrium of Judd&rsquo;s model can be computed by filling in the matrices appropriately.\n",
    "\n",
    "The exercise is to calculate these matrices and compute the following figures.\n",
    "\n",
    "The first figure shows the dynamics of inventories for each firm when the parameters are"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "&delta; = 0.02\n",
    "D = np.array([[-1, 0.5], [0.5, -1]])\n",
    "b = np.array([25, 25])\n",
    "c1 = c2 = np.array([1, -2, 1])\n",
    "e1 = e2 = np.array([10, 10, 3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://s3-ap-southeast-2.amazonaws.com/lectures.quantecon.org/py/_static/lecture_specific/markov_perf/judd_fig2.png\" style=\"\">\n",
    "\n",
    "  \n",
    "Inventories trend to a common steady state.\n",
    "\n",
    "If we increase the depreciation rate to $ \\delta = 0.05 $, then we expect steady state inventories to fall.\n",
    "\n",
    "This is indeed the case, as the next figure shows\n",
    "\n",
    "<img src=\"https://s3-ap-southeast-2.amazonaws.com/lectures.quantecon.org/py/_static/lecture_specific/markov_perf/judd_fig1.png\" style=\"\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "First, let&rsquo;s compute the duopoly MPE under the stated parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# == Parameters == #\n",
    "a0 = 10.0\n",
    "a1 = 2.0\n",
    "&beta; = 0.96\n",
    "&gamma; = 12.0\n",
    "\n",
    "# == In LQ form == #\n",
    "A  = np.eye(3)\n",
    "B1 = np.array([[0.], [1.], [0.]])\n",
    "B2 = np.array([[0.], [0.], [1.]])\n",
    "R1 = [[      0.,      -a0/2,          0.],\n",
    "      [-a0 / 2.,         a1,     a1 / 2.],\n",
    "      [       0,    a1 / 2.,          0.]]\n",
    "\n",
    "R2 = [[     0.,          0.,     -a0 / 2],\n",
    "      [     0.,          0.,     a1 / 2.],\n",
    "      [-a0 / 2,     a1 / 2.,          a1]]\n",
    "\n",
    "Q1 = Q2 = &gamma;\n",
    "S1 = S2 = W1 = W2 = M1 = M2 = 0.0\n",
    "\n",
    "# == Solve using QE's nnash function == #\n",
    "F1, F2, P1, P2 = qe.nnash(A, B1, B2, R1, R2, Q1,\n",
    "                          Q2, S1, S2, W1, W2, M1,\n",
    "                          M2, beta=&beta;)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we evaluate the time path of industry output and prices given\n",
    "initial condition $ q_{10} = q_{20} = 1 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "AF = A - B1 @ F1 - B2 @ F2\n",
    "n = 20\n",
    "x = np.empty((3, n))\n",
    "x[:, 0] = 1, 1, 1\n",
    "for t in range(n-1):\n",
    "    x[:, t+1] = AF @ x[:, t]\n",
    "q1 = x[1, :]\n",
    "q2 = x[2, :]\n",
    "q = q1 + q2       # Total output, MPE\n",
    "p = a0 - a1 * q   # Price, MPE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let&rsquo;s have a look at the monopoly solution.\n",
    "\n",
    "For the state and control, we take\n",
    "\n",
    "$$\n",
    "x_t = q_t - \\bar q\n",
    "\\quad \\text{and} \\quad\n",
    "u_t = q_{t+1} - q_t\n",
    "$$\n",
    "\n",
    "To convert to an LQ problem we set\n",
    "\n",
    "$$\n",
    "R = a_1\n",
    "\\quad \\text{and} \\quad\n",
    "Q = \\gamma\n",
    "$$\n",
    "\n",
    "in the payoff function $ x_t' R x_t + u_t' Q u_t $ and\n",
    "\n",
    "$$\n",
    "A = B = 1\n",
    "$$\n",
    "\n",
    "in the law of motion $ x_{t+1} = A x_t + B u_t $.\n",
    "\n",
    "We solve for the optimal policy $ u_t = - Fx_t $ and track the\n",
    "resulting dynamics of $ \\{q_t\\} $, starting at $ q_0 = 2.0 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "R = a1\n",
    "Q = &gamma;\n",
    "A = B = 1\n",
    "lq_alt = qe.LQ(Q, R, A, B, beta=&beta;)\n",
    "P, F, d = lq_alt.stationary_values()\n",
    "q_bar = a0 / (2.0 * a1)\n",
    "qm = np.empty(n)\n",
    "qm[0] = 2\n",
    "x0 = qm[0] - q_bar\n",
    "x = x0\n",
    "for i in range(1, n):\n",
    "    x = A * x - B * F * x\n",
    "    qm[i] = float(x) + q_bar\n",
    "pm = a0 - a1 * qm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let&rsquo;s have a look at the different time paths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 1, figsize=(9, 9))\n",
    "\n",
    "ax = axes[0]\n",
    "ax.plot(qm, 'b-', lw=2, alpha=0.75, label='monopolist output')\n",
    "ax.plot(q, 'g-', lw=2, alpha=0.75, label='MPE total output')\n",
    "ax.set(ylabel=\"output\", xlabel=\"time\", ylim=(2, 4))\n",
    "ax.legend(loc='upper left', frameon=0)\n",
    "\n",
    "ax = axes[1]\n",
    "ax.plot(pm, 'b-', lw=2, alpha=0.75, label='monopolist price')\n",
    "ax.plot(p, 'g-', lw=2, alpha=0.75, label='MPE price')\n",
    "ax.set(ylabel=\"price\", xlabel=\"time\")\n",
    "ax.legend(loc='upper right', frameon=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2\n",
    "\n",
    "We treat the case $ \\delta = 0.02 $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "&delta; = 0.02\n",
    "D = np.array([[-1, 0.5], [0.5, -1]])\n",
    "b = np.array([25, 25])\n",
    "c1 = c2 = np.array([1, -2, 1])\n",
    "e1 = e2 = np.array([10, 10, 3])\n",
    "\n",
    "&delta;_1 = 1 - &delta;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recalling that the control and state are\n",
    "\n",
    "$$\n",
    "u_{it} =\n",
    "\\begin{bmatrix}\n",
    "    p_{it} \\\\\n",
    "    q_{it}\n",
    "\\end{bmatrix}\n",
    "\\quad \\text{and} \\quad\n",
    "x_t =\n",
    "\\begin{bmatrix}\n",
    "    I_{1t} \\\\\n",
    "    I_{2t} \\\\\n",
    "    1\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "we set up the matrices as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# ==  Create matrices needed to compute the Nash feedback equilibrium == #\n",
    "\n",
    "A = np.array([[&delta;_1,      0,    -&delta;_1 * b[0]],\n",
    "              [  0,    &delta;_1,    -&delta;_1 * b[1]],\n",
    "              [  0,      0,             1]])\n",
    "\n",
    "B1 = &delta;_1 * np.array([[1, -D[0, 0]],\n",
    "                     [0, -D[1, 0]],\n",
    "                     [0,       0]])\n",
    "B2 = &delta;_1 * np.array([[0, -D[0, 1]],\n",
    "                     [1, -D[1, 1]],\n",
    "                     [0,       0]])\n",
    "\n",
    "R1 = -np.array([[0.5 * c1[2],     0,    0.5 * c1[1]],\n",
    "                [          0,     0,              0],\n",
    "                [0.5 * c1[1],     0,         c1[0]]])\n",
    "R2 = -np.array([[0,               0,              0],\n",
    "                [0,     0.5 * c2[2],    0.5 * c2[1]],\n",
    "                [0,     0.5 * c2[1],         c2[0]]])\n",
    "\n",
    "Q1 = np.array([[-0.5 * e1[2], 0], [0, D[0, 0]]])\n",
    "Q2 = np.array([[-0.5 * e2[2], 0], [0, D[1, 1]]])\n",
    "\n",
    "S1 = np.zeros((2, 2))\n",
    "S2 = np.copy(S1)\n",
    "\n",
    "W1 = np.array([[           0,             0],\n",
    "               [           0,             0],\n",
    "               [-0.5 * e1[1],    b[0] / 2.]])\n",
    "W2 = np.array([[           0,             0],\n",
    "               [           0,             0],\n",
    "               [-0.5 * e2[1],    b[1] / 2.]])\n",
    "\n",
    "M1 = np.array([[0, 0], [0, D[0, 1] / 2.]])\n",
    "M2 = np.copy(M1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now compute the equilibrium using `qe.nnash`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "F1, F2, P1, P2 = qe.nnash(A, B1, B2, R1,\n",
    "                          R2, Q1, Q2, S1,\n",
    "                          S2, W1, W2, M1, M2)\n",
    "\n",
    "print(\"\\nFirm 1's feedback rule:\\n\")\n",
    "print(F1)\n",
    "\n",
    "print(\"\\nFirm 2's feedback rule:\\n\")\n",
    "print(F2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let&rsquo;s look at the dynamics of inventories, and reproduce the graph\n",
    "corresponding to $ \\delta = 0.02 $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "AF = A - B1 @ F1 - B2 @ F2\n",
    "n = 25\n",
    "x = np.empty((3, n))\n",
    "x[:, 0] = 2, 0, 1\n",
    "for t in range(n-1):\n",
    "    x[:, t+1] = AF @ x[:, t]\n",
    "I1 = x[0, :]\n",
    "I2 = x[1, :]\n",
    "fig, ax = plt.subplots(figsize=(9, 5))\n",
    "ax.plot(I1, 'b-', lw=2, alpha=0.75, label='inventories, firm 1')\n",
    "ax.plot(I2, 'g-', lw=2, alpha=0.75, label='inventories, firm 2')\n",
    "ax.set_title(rf'$\\delta = {&delta;}$')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "filename": "markov_perf.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Markov Perfect Equilibrium"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}