rob_markov_perf.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e0fae103",
   "metadata": {},
   "source": [
    "\n",
    "<a id='rob-markov-perf'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a080ee65",
   "metadata": {},
   "source": [
    "# Robust Markov Perfect Equilibrium\n",
    "\n",
    "In addition to what’s in Anaconda, this lecture will need the following libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d69bda0b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "!pip install --upgrade quantecon"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21aae7d0",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "\n",
    "<a id='index-0'></a>\n",
    "This lecture describes  a Markov perfect equilibrium  with robust agents.\n",
    "\n",
    "We  focus on special settings with\n",
    "\n",
    "- two players  \n",
    "- quadratic payoff functions  \n",
    "- linear transition rules for the state vector  \n",
    "\n",
    "\n",
    "These specifications simplify calculations and allow us to give a simple example that illustrates basic forces.\n",
    "\n",
    "This lecture is based on ideas described in chapter 15 of  [[Hansen and Sargent, 2008](https://python-advanced.quantecon.org/zreferences.html#id158)] and in [Markov perfect equilibrium](https://python-intro.quantecon.org/markov_perf.html)\n",
    "and [Robustness](https://python-advanced.quantecon.org/robustness.html).\n",
    "\n",
    "Let’s start with some standard imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "767ad13f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import quantecon as qe\n",
    "from scipy.linalg import solve\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6c53253",
   "metadata": {},
   "source": [
    "### Basic Setup\n",
    "\n",
    "Decisions of two agents affect the motion of a  state vector\n",
    "that appears as an argument of   payoff functions of both agents.\n",
    "\n",
    "As described in [Markov perfect equilibrium](https://python-intro.quantecon.org/markov_perf.html), when decision-makers have no concerns about the robustness of their decision rules\n",
    "to misspecifications of the state dynamics, a Markov perfect equilibrium  can be computed via\n",
    "backward recursion  on  two sets of equations\n",
    "\n",
    "- a pair of Bellman equations, one for each agent.  \n",
    "- a pair of equations that express linear decision rules for each agent as functions of that agent’s  continuation value function as well as parameters of preferences and state transition matrices.  \n",
    "\n",
    "\n",
    "This lecture shows how  a similar equilibrium concept and similar computational  procedures\n",
    "apply when we impute concerns about robustness to both decision-makers.\n",
    "\n",
    "A Markov perfect equilibrium with robust agents will be characterized by\n",
    "\n",
    "- a pair of Bellman equations, one for each agent.  \n",
    "- a pair of equations that express linear decision rules for each agent as functions of that agent’s  continuation value function as well as parameters of\n",
    "  preferences and state transition matrices.  \n",
    "- a pair of equations that express linear decision rules for worst-case shocks for each agent as functions of that agent’s  continuation value function as well as parameters of\n",
    "  preferences and state transition matrices.  \n",
    "\n",
    "\n",
    "Below, we’ll construct a robust firms version of the classic duopoly model with\n",
    "adjustment costs analyzed in [Markov perfect equilibrium](https://python-intro.quantecon.org/markov_perf.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3dd57ab",
   "metadata": {},
   "source": [
    "## Linear Markov Perfect Equilibria with Robust Agents\n",
    "\n",
    "\n",
    "<a id='index-1'></a>\n",
    "As we saw in [Markov perfect equilibrium](https://python-intro.quantecon.org/markov_perf.html), the study of Markov perfect equilibria in dynamic games with two players\n",
    "leads us to an interrelated pair of Bellman equations.\n",
    "\n",
    "In linear quadratic dynamic games, these “stacked Bellman equations” become “stacked Riccati equations” with a tractable mathematical structure."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d66f30be",
   "metadata": {},
   "source": [
    "### Modified Coupled Linear Regulator Problems\n",
    "\n",
    "We consider a general linear quadratic regulator game with two players, each of whom fears model misspecifications.\n",
    "\n",
    "We often call the players agents.\n",
    "\n",
    "The agents share a common baseline model for the transition dynamics of the state vector\n",
    "\n",
    "- this is a counterpart of a ‘rational expectations’ assumption of shared beliefs  \n",
    "\n",
    "\n",
    "But now one or more agents doubt that the baseline model is correctly specified.\n",
    "\n",
    "The agents express the possibility that their baseline specification is incorrect by adding a contribution $ C v_{it} $ to the time $ t $ transition law for the state\n",
    "\n",
    "- $ C $ is the usual *volatility matrix* that appears in stochastic versions of optimal linear regulator problems.  \n",
    "- $ v_{it} $ is a possibly history-dependent vector of distortions to the dynamics of the state that agent $ i $ uses to  represent misspecification of the original model.  \n",
    "\n",
    "\n",
    "For convenience, we’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 a sequence $ \\{u_{-it}\\} $ as given and chooses a sequence $ \\{u_{it}\\} $ to minimize and $ \\{v_{it}\\} $ to maximize\n",
    "\n",
    "\n",
    "<a id='equation-rmp-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",
    "    \\theta_i v_{it}' v_{it}\n",
    "\\right\\} \\tag{27.1}\n",
    "$$\n",
    "\n",
    "while  thinking that the state evolves according to\n",
    "\n",
    "\n",
    "<a id='equation-rmp-orig-0'></a>\n",
    "$$\n",
    "x_{t+1} = A x_t + B_1 u_{1t} + B_2 u_{2t} + C v_{it} \\tag{27.2}\n",
    "$$\n",
    "\n",
    "Here\n",
    "\n",
    "- $ x_t $ is an $ n \\times 1 $ state vector,   $ u_{it} $ is a $ k_i \\times 1 $ vector of controls for player $ i $, and  \n",
    "- $ v_{it} $ is an $ h \\times 1 $ vector of distortions to the state dynamics that concern  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 $  \n",
    "- $ C $ is $ n \\times h $  \n",
    "- $ \\theta_i \\in [\\underline \\theta_i, +\\infty] $ is a scalar multiplier parameter of player $ i $  \n",
    "\n",
    "\n",
    "If $ \\theta_i = + \\infty $, player $ i $ completely trusts the baseline model.\n",
    "\n",
    "If $ \\theta_i < _\\infty $, player $ i $ suspects that some other unspecified model actually governs the transition dynamics.\n",
    "\n",
    "The term $ \\theta_i v_{it}' v_{it} $ is a time $ t $ contribution to an entropy penalty that an (imaginary) loss-maximizing agent inside\n",
    "agent $ i $’s mind charges for distorting the law of motion in a way that harms agent $ i $.\n",
    "\n",
    "- the imaginary loss-maximizing  agent helps the loss-minimizing agent by helping him construct bounds on the behavior of his decision rule over a\n",
    "  large **set** of alternative models of state transition dynamics.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cab273c9",
   "metadata": {},
   "source": [
    "### Computing Equilibrium\n",
    "\n",
    "We formulate a linear robust 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",
    "Player $ i $’s malevolent alter ego employs decision rules $ v_{it} = K_{it} x_t $ where $ K_{it} $ is an  $ h \\times n $ matrix.\n",
    "\n",
    "A robust Markov perfect equilibrium is a pair of sequences $ \\{F_{1t}, F_{2t}\\} $ and  a pair of sequences $ \\{K_{1t}, K_{2t}\\} $ over $ t = t_0, \\ldots, t_1 - 1 $ that satisfy\n",
    "\n",
    "- $ \\{F_{1t}, K_{1t}\\} $ solves player 1’s robust decision problem, taking $ \\{F_{2t}\\} $ as given, and  \n",
    "- $ \\{F_{2t}, K_{2t}\\} $ solves player 2’s robust decision problem, taking $ \\{F_{1t}\\} $ as given.  \n",
    "\n",
    "\n",
    "If we substitute $ u_{2t} = - F_{2t} x_t $  into [(27.1)](#equation-rmp-orig-1) and [(27.2)](#equation-rmp-orig-0), then player 1’s problem becomes minimization-maximization of\n",
    "\n",
    "\n",
    "<a id='equation-rmp-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",
    "    \\theta_1 v_{1t}' v_{1t}\n",
    "    \\right\\} \\tag{27.3}\n",
    "$$\n",
    "\n",
    "subject to\n",
    "\n",
    "\n",
    "<a id='equation-rmp-eq-mpe-p1d'></a>\n",
    "$$\n",
    "x_{t+1} = \\Lambda_{1t} x_t + B_1 u_{1t} + C v_{1t} \\tag{27.4}\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 robust dynamic programming problem of the type studied in the [Robustness](https://python-advanced.quantecon.org/robustness.html) lecture,\n",
    "which can be solved by working backward.\n",
    "\n",
    "Maximization with respect to distortion  $ v_{1t} $ leads to  the following version of the $ \\mathcal D $ operator from\n",
    "the [Robustness](https://python-advanced.quantecon.org/robustness.html) lecture, namely\n",
    "\n",
    "\n",
    "<a id='equation-rmp-rb-d'></a>\n",
    "$$\n",
    "\\mathcal D_1(P) := P + PC (\\theta_1 I - C' P C)^{-1} C' P \\tag{27.5}\n",
    "$$\n",
    "\n",
    "The matrix $ F_{1t} $ in the policy rule $ u_{1t} = - F_{1t} x_t $ that solves agent 1’s  problem satisfies\n",
    "\n",
    "\n",
    "<a id='equation-rmp-orig-3'></a>\n",
    "$$\n",
    "F_{1t}\n",
    "= (Q_1 + \\beta B_1' {\\mathcal D}_1( P_{1t+1}) B_1)^{-1}\n",
    "(\\beta B_1' {\\mathcal D}_1( P_{1t+1}) \\Lambda_{1t} + \\Gamma_{1t}) \\tag{27.6}\n",
    "$$\n",
    "\n",
    "where $ P_{1t} $ solves the matrix Riccati difference equation\n",
    "\n",
    "\n",
    "<a id='equation-rmp-orig-4'></a>\n",
    "$$\n",
    "P_{1t} =\n",
    "\\Pi_{1t} -\n",
    "(\\beta B_1' {\\mathcal D}_1(P_{1t+1}) \\Lambda_{1t} + \\Gamma_{1t})' (Q_1 + \\beta B_1' {\\mathcal D}_1( P_{1t+1}) B_1)^{-1}\n",
    "(\\beta B_1' {\\mathcal D}_1(P_{1t+1}) \\Lambda_{1t} + \\Gamma_{1t}) + \\\\\n",
    "\\beta \\Lambda_{1t}' {\\mathcal D}_1(P_{1t+1}) \\Lambda_{1t} \\tag{27.7}\n",
    "$$\n",
    "\n",
    "Similarly, the policy that solves player 2’s problem is\n",
    "\n",
    "\n",
    "<a id='equation-rmp-orig-5'></a>\n",
    "$$\n",
    "F_{2t} = (Q_2 + \\beta B_2' {\\mathcal D}_2( P_{2t+1} ) B_2)^{-1}\n",
    "(\\beta B_2' {\\mathcal D}_2 (P_{2t+1}) \\Lambda_{2t} + \\Gamma_{2t}) \\tag{27.8}\n",
    "$$\n",
    "\n",
    "where $ P_{2t} $ solves\n",
    "\n",
    "\n",
    "<a id='equation-rmp-orig-6'></a>\n",
    "$$\n",
    "P_{2t} =\n",
    "\\Pi_{2t} - (\\beta B_2' {\\mathcal D}_2 ( P_{2t+1}) \\Lambda_{2t} + \\Gamma_{2t})' (Q_2 + \\beta B_2' {\\mathcal D}_2 ( P_{2t+1}) B_2)^{-1}(\\beta B_2' {\\mathcal D}_2 ( P_{2t+1}) \\Lambda_{2t} + \\Gamma_{2t}) + \\\\\n",
    "\\beta \\Lambda_{2t}' {\\mathcal D}_2 ( P_{2t+1})  \\Lambda_{2t} \\tag{27.9}\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 [(27.6)](#equation-rmp-orig-3), [(27.7)](#equation-rmp-orig-4), [(27.8)](#equation-rmp-orig-5), and [(27.9)](#equation-rmp-orig-6), and “work backwards” from time $ t_1 - 1 $.\n",
    "\n",
    "Since we’re working backwards, $ 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 [(27.6)](#equation-rmp-orig-3) contain $ F_{2t} $  \n",
    "- some terms on the right-hand side of [(27.8)](#equation-rmp-orig-5) contain $ F_{1t} $  \n",
    "\n",
    "\n",
    "we need to solve these $ k_1 + k_2 $ equations simultaneously."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "253f7a11",
   "metadata": {},
   "source": [
    "### Key Insight\n",
    "\n",
    "As in [Markov perfect equilibrium](https://python-intro.quantecon.org/markov_perf.html), a key insight here  is that  equations  [(27.6)](#equation-rmp-orig-3) and [(27.8)](#equation-rmp-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 [(27.7)](#equation-rmp-orig-4) and [(27.9)](#equation-rmp-orig-6).\n",
    "\n",
    "Notice how $ j $’s control law $ F_{jt} $ is a function of $ \\{F_{is}, s \\geq t, i \\neq j \\} $.\n",
    "\n",
    "Thus, agent $ i $’s choice of $ \\{F_{it}; t = t_0, \\ldots, t_1 - 1\\} $ influences agent $ j $’s choice of control laws.\n",
    "\n",
    "However, in the Markov perfect equilibrium of this game, each agent is assumed to ignore the influence that his choice exerts on the other agent’s choice.\n",
    "\n",
    "After these equations have been solved, we can also deduce associated sequences of worst-case shocks."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ec1592a",
   "metadata": {},
   "source": [
    "### Worst-case Shocks\n",
    "\n",
    "For agent $ i $ the maximizing or worst-case shock $ v_{it} $ is\n",
    "\n",
    "$$\n",
    "v_{it} = K_{it} x_t\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "K_{it} = \\theta_i^{-1} (I - \\theta_i^{-1} C' P_{i,t+1} C)^{-1} C' P_{i,t+1} (A - B_1 F_{it} - B_2 F_{2t})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a94fc868",
   "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",
   "id": "c6af9d5d",
   "metadata": {},
   "source": [
    "### Implementation\n",
    "\n",
    "We use the function nnash_robust to compute a\n",
    "Markov perfect equilibrium of the infinite horizon linear quadratic dynamic\n",
    "game with robust planers in the manner described above."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39c460bf",
   "metadata": {},
   "source": [
    "## Application\n",
    "\n",
    "\n",
    "<a id='index-2'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d5f332b",
   "metadata": {},
   "source": [
    "### A Duopoly Model\n",
    "\n",
    "Without concerns for robustness, the model is identical to the duopoly model from the [Markov perfect equilibrium](https://python-intro.quantecon.org/markov_perf.html) lecture.\n",
    "\n",
    "To begin, we briefly review the structure of that 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-rgame2'></a>\n",
    "$$\n",
    "p = a_0 - a_1 (q_1 +  q_2) \\tag{27.10}\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 [(27.10)](#equation-rgame2) 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-rgame1'></a>\n",
    "$$\n",
    "\\pi_i = p q_i - \\gamma (\\hat q_i - q_i)^2, \\quad \\gamma > 0 , \\tag{27.11}\n",
    "$$\n",
    "\n",
    "Substituting the inverse demand curve [(27.10)](#equation-rgame2) into [(27.11)](#equation-rgame1) lets us express the one-period payoff as\n",
    "\n",
    "\n",
    "<a id='equation-rgame3'></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{27.12}\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",
    "This completes our review of the duopoly model without concerns for robustness.\n",
    "\n",
    "Now we activate robustness  concerns of both firms.\n",
    "\n",
    "To map a robust version of the duopoly model  into  coupled robust linear-quadratic dynamic programming\n",
    "problems, we again define the state 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 [(27.11)](#equation-rgame1) for the two firms in the duopoly model.\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",
    "A robust  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-rmp-eq-mpe-cle'></a>\n",
    "$$\n",
    "x_{t+1} = (A - B_1 F_1 -B_1 F_2 ) x_t \\tag{27.13}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c058158a",
   "metadata": {},
   "source": [
    "### Parameters and Solution\n",
    "\n",
    "Consider the  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 computed the infinite horizon MPE without robustness using the code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31776d65",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import quantecon as qe\n",
    "\n",
    "# Parameters\n",
    "a0 = 10.0\n",
    "a1 = 2.0\n",
    "β = 0.96\n",
    "γ = 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 = γ\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=β)\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",
   "id": "bc043210",
   "metadata": {},
   "source": [
    "#### Markov Perfect Equilibrium with Robustness\n",
    "\n",
    "We add robustness concerns to the Markov Perfect Equilibrium model by\n",
    "extending the function `qe.nnash`\n",
    "([link](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/lqnash.py))\n",
    "into a robustness version by adding the maximization operator\n",
    "$ \\mathcal D(P) $ into the backward induction.\n",
    "\n",
    "The MPE with robustness function is `nnash_robust`.\n",
    "\n",
    "The function’s code is as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "264f31ce",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def nnash_robust(A, C, B1, B2, R1, R2, Q1, Q2, S1, S2, W1, W2, M1, M2,\n",
    "                 θ1, θ2, beta=1.0, tol=1e-8, max_iter=1000):\n",
    "\n",
    "    \"\"\"\n",
    "    Compute the limit of a Nash linear quadratic dynamic game with\n",
    "    robustness concern.\n",
    "\n",
    "    In this problem, player i minimizes\n",
    "    .. math::\n",
    "        \\sum_{t=0}^{\\infty}\n",
    "        \\left\\{\n",
    "            x_t' r_i x_t + 2 x_t' w_i\n",
    "            u_{it} +u_{it}' q_i u_{it} + u_{jt}' s_i u_{jt} + 2 u_{jt}'\n",
    "            m_i u_{it}\n",
    "        \\right\\}\n",
    "    subject to the law of motion\n",
    "    .. math::\n",
    "        x_{it+1} = A x_t + b_1 u_{1t} + b_2 u_{2t} + C w_{it+1}\n",
    "    and a perceived control law :math:`u_j(t) = - f_j x_t` for the other\n",
    "    player.\n",
    "\n",
    "    The player i also concerns about the model misspecification,\n",
    "    and maximizes\n",
    "    .. math::\n",
    "        \\sum_{t=0}^{\\infty}\n",
    "        \\left\\{\n",
    "            \\beta^{t+1} \\theta_{i} w_{it+1}'w_{it+1}\n",
    "        \\right\\}\n",
    "\n",
    "    The solution computed in this routine is the :math:`f_i` and\n",
    "    :math:`P_i` of the associated double optimal linear regulator\n",
    "    problem.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    A : scalar(float) or array_like(float)\n",
    "        Corresponds to the MPE equations, should be of size (n, n)\n",
    "    C : scalar(float) or array_like(float)\n",
    "        As above, size (n, c), c is the size of w\n",
    "    B1 : scalar(float) or array_like(float)\n",
    "        As above, size (n, k_1)\n",
    "    B2 : scalar(float) or array_like(float)\n",
    "        As above, size (n, k_2)\n",
    "    R1 : scalar(float) or array_like(float)\n",
    "        As above, size (n, n)\n",
    "    R2 : scalar(float) or array_like(float)\n",
    "        As above, size (n, n)\n",
    "    Q1 : scalar(float) or array_like(float)\n",
    "        As above, size (k_1, k_1)\n",
    "    Q2 : scalar(float) or array_like(float)\n",
    "        As above, size (k_2, k_2)\n",
    "    S1 : scalar(float) or array_like(float)\n",
    "        As above, size (k_1, k_1)\n",
    "    S2 : scalar(float) or array_like(float)\n",
    "        As above, size (k_2, k_2)\n",
    "    W1 : scalar(float) or array_like(float)\n",
    "        As above, size (n, k_1)\n",
    "    W2 : scalar(float) or array_like(float)\n",
    "        As above, size (n, k_2)\n",
    "    M1 : scalar(float) or array_like(float)\n",
    "        As above, size (k_2, k_1)\n",
    "    M2 : scalar(float) or array_like(float)\n",
    "        As above, size (k_1, k_2)\n",
    "    θ1 : scalar(float)\n",
    "            Robustness parameter of player 1\n",
    "    θ2 : scalar(float)\n",
    "            Robustness parameter of player 2\n",
    "    beta : scalar(float), optional(default=1.0)\n",
    "          Discount factor\n",
    "    tol : scalar(float), optional(default=1e-8)\n",
    "         This is the tolerance level for convergence\n",
    "    max_iter : scalar(int), optional(default=1000)\n",
    "        This is the maximum number of iterations allowed\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    F1 : array_like, dtype=float, shape=(k_1, n)\n",
    "        Feedback law for agent 1\n",
    "    F2 : array_like, dtype=float, shape=(k_2, n)\n",
    "        Feedback law for agent 2\n",
    "    P1 : array_like, dtype=float, shape=(n, n)\n",
    "        The steady-state solution to the associated discrete matrix\n",
    "        Riccati equation for agent 1\n",
    "    P2 : array_like, dtype=float, shape=(n, n)\n",
    "        The steady-state solution to the associated discrete matrix\n",
    "        Riccati equation for agent 2\n",
    "    \"\"\"\n",
    "\n",
    "    # Unload parameters and make sure everything is a matrix\n",
    "    params = A, C, B1, B2, R1, R2, Q1, Q2, S1, S2, W1, W2, M1, M2\n",
    "    params = map(np.asmatrix, params)\n",
    "    A, C, B1, B2, R1, R2, Q1, Q2, S1, S2, W1, W2, M1, M2 = params\n",
    "\n",
    "\n",
    "    # Multiply A, B1, B2 by sqrt(β) to enforce discounting\n",
    "    A, B1, B2 = [np.sqrt(β) * x for x in (A, B1, B2)]\n",
    "\n",
    "    # Initial values\n",
    "    n = A.shape[0]\n",
    "    k_1 = B1.shape[1]\n",
    "    k_2 = B2.shape[1]\n",
    "\n",
    "    v1 = np.eye(k_1)\n",
    "    v2 = np.eye(k_2)\n",
    "    P1 = np.eye(n) * 1e-5\n",
    "    P2 = np.eye(n) * 1e-5\n",
    "    F1 = np.random.randn(k_1, n)\n",
    "    F2 = np.random.randn(k_2, n)\n",
    "\n",
    "\n",
    "    for it in range(max_iter):\n",
    "        # Update\n",
    "        F10 = F1\n",
    "        F20 = F2\n",
    "\n",
    "        I = np.eye(C.shape[1])\n",
    "\n",
    "        # D1(P1)\n",
    "        # Note: INV1 may not be solved if the matrix is singular\n",
    "        INV1 = solve(θ1 * I - C.T @ P1 @ C, I)\n",
    "        D1P1 =  P1 + P1 @ C @ INV1 @ C.T @ P1\n",
    "\n",
    "        # D2(P2)\n",
    "        # Note: INV2 may not be solved if the matrix is singular\n",
    "        INV2 = solve(θ2 * I - C.T @ P2 @ C, I)\n",
    "        D2P2 =  P2 + P2 @ C @ INV2 @ C.T @ P2\n",
    "\n",
    "        G2 = solve(Q2 + B2.T @ D2P2 @ B2, v2)\n",
    "        G1 = solve(Q1 + B1.T @ D1P1 @ B1, v1)\n",
    "        H2 = G2 @ B2.T @ D2P2\n",
    "        H1 = G1 @ B1.T @ D1P1\n",
    "\n",
    "        # Break up the computation of F1, F2\n",
    "        F1_left = v1 - (H1 @ B2 + G1 @ M1.T) @ (H2 @ B1 + G2 @ M2.T)\n",
    "        F1_right = H1 @ A + G1 @ W1.T - \\\n",
    "                   (H1 @ B2 + G1 @ M1.T) @ (H2 @ A + G2 @ W2.T)\n",
    "        F1 = solve(F1_left, F1_right)\n",
    "        F2 = H2 @ A + G2 @ W2.T - (H2 @ B1 + G2 @ M2.T) @ F1\n",
    "\n",
    "        Λ1 = A - B2 @ F2\n",
    "        Λ2 = A - B1 @ F1\n",
    "        Π1 = R1 + F2.T @ S1 @ F2\n",
    "        Π2 = R2 + F1.T @ S2 @ F1\n",
    "        Γ1 = W1.T - M1.T @ F2\n",
    "        Γ2 = W2.T - M2.T @ F1\n",
    "\n",
    "        # Compute P1 and P2\n",
    "        P1 = Π1 - (B1.T @ D1P1 @ Λ1 + Γ1).T @ F1 + \\\n",
    "             Λ1.T @ D1P1 @ Λ1\n",
    "        P2 = Π2 - (B2.T @ D2P2 @ Λ2 + Γ2).T @ F2 + \\\n",
    "             Λ2.T @ D2P2 @ Λ2\n",
    "\n",
    "        dd = np.max(np.abs(F10 - F1)) + np.max(np.abs(F20 - F2))\n",
    "\n",
    "        if dd < tol:  # success!\n",
    "            break\n",
    "\n",
    "    else:\n",
    "        raise ValueError(f'No convergence: Iteration limit of {max_iter} \\\n",
    "            reached in nnash')\n",
    "\n",
    "    return F1, F2, P1, P2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a39d480",
   "metadata": {},
   "source": [
    "### Some Details\n",
    "\n",
    "Firm  $ i $ wants to minimize\n",
    "\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\\}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\begin{split}\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",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "\\begin{split}\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",
    "\\text{,} \\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",
    "\\text{,} \\quad\n",
    "Q_1 = Q_2 = \\gamma\n",
    "\\text{,} \\quad\n",
    "S_1 = S_2 = 0\n",
    "\\text{,} \\quad\n",
    "W_1 = W_2 = 0\n",
    "\\text{,} \\quad\n",
    "M_1 = M_2 = 0\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "The parameters of the duopoly model are:\n",
    "\n",
    "- $ a_0 = 10 $  \n",
    "- $ a_1 = 2 $  \n",
    "- $ \\beta = 0.96 $  \n",
    "- $ \\gamma = 12 $  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64416b72",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Parameters\n",
    "a0 = 10.0\n",
    "a1 = 2.0\n",
    "β = 0.96\n",
    "γ = 12.0\n",
    "\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 = γ\n",
    "S1 = S2 = W1 = W2 = M1 = M2 = 0.0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0e47961",
   "metadata": {},
   "source": [
    "#### Consistency Check\n",
    "\n",
    "We first conduct a comparison test to check if `nnash_robust` agrees\n",
    "with `qe.nnash` in the non-robustness case in which each $ \\theta_i \\approx +\\infty $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91722bed",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# 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=β)\n",
    "\n",
    "# Solve using nnash_robust\n",
    "F1r, F2r, P1r, P2r = nnash_robust(A, np.zeros((3, 1)), B1, B2, R1, R2, Q1,\n",
    "                                  Q2, S1, S2, W1, W2, M1, M2, 1e-10,\n",
    "                                  1e-10, beta=β)\n",
    "\n",
    "print('F1 and F1r should be the same: ', np.allclose(F1, F1r))\n",
    "print('F2 and F2r should be the same: ', np.allclose(F1, F1r))\n",
    "print('P1 and P1r should be the same: ', np.allclose(P1, P1r))\n",
    "print('P2 and P2r should be the same: ', np.allclose(P1, P1r))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b37b0e3",
   "metadata": {},
   "source": [
    "We can see that the results are consistent across the two functions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dbbb854b",
   "metadata": {},
   "source": [
    "#### Comparative Dynamics under Baseline Transition Dynamics\n",
    "\n",
    "We want to compare the dynamics of price and output under the baseline\n",
    "MPE model with those under the baseline model under the robust decision rules within the robust MPE.\n",
    "\n",
    "This means that we simulate the state dynamics under the MPE  equilibrium **closed-loop** transition matrix\n",
    "\n",
    "$$\n",
    "A^o = A - B_1 F_1 - B_2 F_2\n",
    "$$\n",
    "\n",
    "where $ F_1 $ and $ F_2 $ are the firms’ robust decision rules within the robust markov_perfect equilibrium\n",
    "\n",
    "- by simulating under the baseline model transition dynamics and the robust MPE rules we are in assuming that at the end of the day\n",
    "  firms’ concerns about misspecification of the baseline model do not materialize.  \n",
    "- a short way of saying this is that misspecification fears are all ‘just in the minds’ of the firms.  \n",
    "- simulating under the baseline model is a common practice in the literature.  \n",
    "- note that *some* assumption about the model that actually governs the data has to be made in order to create a simulation.  \n",
    "- later we will describe the (erroneous) beliefs of the two firms  that justify their robust decisions as best responses to transition\n",
    "  laws that are distorted relative to the baseline model.  \n",
    "\n",
    "\n",
    "After simulating $ x_t $ under the baseline transition dynamics and robust decision rules $ F_i, i = 1, 2 $, we\n",
    "extract and plot industry output $ q_t=q_{1t}+q_{2t} $ and price $ p_t = a_0 − a_1 q_t $.\n",
    "\n",
    "Here we set the robustness and volatility matrix parameters as follows:\n",
    "\n",
    "- $ \\theta_1 = 0.02 $  \n",
    "- $ \\theta_2 = 0.04 $  \n",
    "- $ C = \\begin{pmatrix} 0 \\\\ 0.01 \\\\ 0.01 \\end{pmatrix} $  \n",
    "\n",
    "\n",
    "Because we have set $ \\theta_1 < \\theta_2 < + \\infty $ we know that\n",
    "\n",
    "- both firms fear that the baseline specification of the state transition dynamics are incorrect.  \n",
    "- firm $ 1 $ fears misspecification more than firm $ 2 $.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4830243d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Robustness parameters and matrix\n",
    "C = np.asmatrix([[0], [0.01], [0.01]])\n",
    "θ1 = 0.02\n",
    "θ2 = 0.04\n",
    "n = 20\n",
    "\n",
    "\n",
    "# Solve using nnash_robust\n",
    "F1r, F2r, P1r, P2r = nnash_robust(A, C, B1, B2, R1, R2, Q1,\n",
    "                                  Q2, S1, S2, W1, W2, M1, M2,\n",
    "                                  θ1, θ2, beta=β)\n",
    "\n",
    "\n",
    "\n",
    "# MPE output and price\n",
    "AF = A - B1 @ F1 - B2 @ F2\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",
    "\n",
    "# RMPE output and price\n",
    "AO = A - B1 @ F1r - B2 @ F2r\n",
    "xr = np.empty((3, n))\n",
    "xr[:, 0] = 1, 1, 1\n",
    "for t in range(n - 1):\n",
    "    xr[:, t+1] = AO @ xr[:, t]\n",
    "qr1 = xr[1, :]\n",
    "qr2 = xr[2, :]\n",
    "qr = qr1 + qr2      # Total output, RMPE\n",
    "pr = a0 - a1 * qr   # Price, RMPE\n",
    "\n",
    "# RMPE heterogeneous beliefs output and price\n",
    "I = np.eye(C.shape[1])\n",
    "INV1 = solve(θ1 * I - C.T @ P1 @ C, I)\n",
    "K1 =  P1 @ C @ INV1 @ C.T @ P1 @ AO\n",
    "AOCK1 = AO + C.T @ K1\n",
    "\n",
    "INV2 = solve(θ2 * I - C.T @ P2 @ C, I)\n",
    "K2 =  P2 @ C @ INV2 @ C.T @ P2 @ AO\n",
    "AOCK2 = AO + C.T @ K2\n",
    "xrp1 = np.empty((3, n))\n",
    "xrp2 = np.empty((3, n))\n",
    "xrp1[:, 0] = 1, 1, 1\n",
    "xrp2[:, 0] = 1, 1, 1\n",
    "for t in range(n - 1):\n",
    "    xrp1[:, t + 1] = AOCK1 @ xrp1[:, t]\n",
    "    xrp2[:, t + 1] = AOCK2 @ xrp2[:, t]\n",
    "qrp11 = xrp1[1, :]\n",
    "qrp12 = xrp1[2, :]\n",
    "qrp21 = xrp2[1, :]\n",
    "qrp22 = xrp2[2, :]\n",
    "qrp1 = qrp11 + qrp12       # Total output, RMPE from player 1's belief\n",
    "qrp2 = qrp21 + qrp22       # Total output, RMPE from player 2's belief\n",
    "prp1 = a0 - a1 * qrp1      # Price, RMPE from player 1's belief\n",
    "prp2 = a0 - a1 * qrp2      # Price, RMPE from player 2's belief"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edfd2570",
   "metadata": {},
   "source": [
    "The following code prepares graphs that compare market-wide output $ q_{1t} + q_{2t} $ and  the price of the good\n",
    "$ p_t $ under equilibrium decision rules $ F_i, i = 1, 2 $ from an ordinary Markov perfect equilibrium and the decision rules\n",
    "under a  Markov perfect equilibrium with robust firms with multiplier parameters $ \\theta_i, i = 1,2 $ set as described above.\n",
    "\n",
    "Both industry output and price  are under the transition dynamics associated with the baseline model; only the decision rules $ F_i $ differ across the two\n",
    "equilibrium objects presented."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c45b80ac",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 1, figsize=(9, 9))\n",
    "\n",
    "ax = axes[0]\n",
    "ax.plot(q, 'g-', lw=2, alpha=0.75, label='MPE output')\n",
    "ax.plot(qr, 'm-', lw=2, alpha=0.75, label='RMPE 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(p, 'g-', lw=2, alpha=0.75, label='MPE price')\n",
    "ax.plot(pr, 'm-', lw=2, alpha=0.75, label='RMPE price')\n",
    "ax.set(ylabel=\"price\", xlabel=\"time\")\n",
    "ax.legend(loc='upper right', frameon=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ec811a1",
   "metadata": {},
   "source": [
    "Under the dynamics associated with the baseline model, the price path is higher with the Markov perfect equilibrium robust decision rules\n",
    "than it is with decision rules for the ordinary Markov perfect equilibrium.\n",
    "\n",
    "So is the industry output path.\n",
    "\n",
    "To dig a little beneath the forces driving these outcomes, we want to plot $ q_{1t} $\n",
    "and $ q_{2t} $ in the Markov perfect equilibrium with robust firms and to compare them with corresponding objects\n",
    "in the Markov perfect equilibrium without robust firms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "675bfede",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 1, figsize=(9, 9))\n",
    "\n",
    "ax = axes[0]\n",
    "ax.plot(q1, 'g-', lw=2, alpha=0.75, label='firm 1 MPE output')\n",
    "ax.plot(qr1, 'b-', lw=2, alpha=0.75, label='firm 1 RMPE output')\n",
    "ax.set(ylabel=\"output\", xlabel=\"time\", ylim=(1, 2))\n",
    "ax.legend(loc='upper left', frameon=0)\n",
    "\n",
    "ax = axes[1]\n",
    "ax.plot(q2, 'g-', lw=2, alpha=0.75, label='firm 2 MPE output')\n",
    "ax.plot(qr2, 'r-', lw=2, alpha=0.75, label='firm 2 RMPE output')\n",
    "ax.set(ylabel=\"output\", xlabel=\"time\", ylim=(1, 2))\n",
    "ax.legend(loc='upper left', frameon=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44fb7c87",
   "metadata": {},
   "source": [
    "Evidently, firm 1’s output path is substantially lower when firms are  robust firms while\n",
    "firm 2’s output path is virtually the same as it would be in an ordinary  Markov perfect equilibrium with no robust firms.\n",
    "\n",
    "Recall that we have set $ \\theta_1 = .02 $ and $ \\theta_2 = .04 $, so that firm 1 fears\n",
    "misspecification of the baseline model substantially more than does firm 2\n",
    "\n",
    "- but also please notice that firm $ 2 $’s  behavior in the Markov perfect equilibrium with robust firms responds to the decision rule $ F_1 x_t $\n",
    "  employed by firm $ 1 $.  \n",
    "- thus it is something of a coincidence that its output is almost the same in the two equilibria.  \n",
    "\n",
    "\n",
    "Larger  concerns about misspecification induce firm 1 to be more cautious than firm 2 in predicting market price and the output of the other firm.\n",
    "\n",
    "To explore this, we study next how *ex-post* the two firms’ beliefs about state dynamics differ in  the Markov perfect equilibrium with robust firms.\n",
    "\n",
    "(by *ex-post* we mean *after* extremization of each firm’s intertemporal objective)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8933a84",
   "metadata": {},
   "source": [
    "#### Heterogeneous Beliefs\n",
    "\n",
    "As before, let $ A^o = A - B\\_1 F\\_1^r - B\\_2 F\\_2^r $, where in a robust MPE, $ F_i^r $ is a robust decision rule for firm $ i $.\n",
    "\n",
    "Worst-case forecasts of  $ x_t $ starting from $ t=0 $ differ between the two firms.\n",
    "\n",
    "This means that  worst-case forecasts of industry output $ q_{1t} + q_{2t} $ and price $ p_t $ also differ between the two firms.\n",
    "\n",
    "To find these worst-case beliefs, we compute the following three “closed-loop” transition matrices\n",
    "\n",
    "- $ A^o $  \n",
    "- $ A^o + C K\\_1 $  \n",
    "- $ A^o + C K\\_2 $  \n",
    "\n",
    "\n",
    "We call the first transition law, namely, $ A^o $,  the baseline transition under firms’ robust decision rules.\n",
    "\n",
    "We call the second and third worst-case transitions under robust decision\n",
    "rules for firms 1 and 2.\n",
    "\n",
    "From $ \\{x_t\\} $ paths generated by  each of these transition laws, we pull off the associated price and total output sequences.\n",
    "\n",
    "The following code plots them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "581d00ac",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "print('Baseline Robust transition matrix AO is: \\n', np.round(AO, 3))\n",
    "print('Player 1\\'s worst-case transition matrix AOCK1 is: \\n', \\\n",
    "np.round(AOCK1, 3))\n",
    "print('Player 2\\'s worst-case transition matrix AOCK2 is: \\n', \\\n",
    "np.round(AOCK2, 3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8b06266",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# == Plot == #\n",
    "fig, axes = plt.subplots(2, 1, figsize=(9, 9))\n",
    "\n",
    "ax = axes[0]\n",
    "ax.plot(qrp1, 'b--', lw=2, alpha=0.75,\n",
    "    label='RMPE worst-case belief output player 1')\n",
    "ax.plot(qrp2, 'r:', lw=2, alpha=0.75,\n",
    "    label='RMPE worst-case belief output player 2')\n",
    "ax.plot(qr, 'm-', lw=2, alpha=0.75, label='RMPE 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(prp1, 'b--', lw=2, alpha=0.75,\n",
    "    label='RMPE worst-case belief price player 1')\n",
    "ax.plot(prp2, 'r:', lw=2, alpha=0.75,\n",
    "    label='RMPE worst-case belief price player 2')\n",
    "ax.plot(pr, 'm-', lw=2, alpha=0.75, label='RMPE price')\n",
    "ax.set(ylabel=\"price\", xlabel=\"time\")\n",
    "ax.legend(loc='upper right', frameon=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e8a6358",
   "metadata": {},
   "source": [
    "We see from the above graph that under robustness concerns, player 1 and\n",
    "player 2  have heterogeneous beliefs about total  output and the goods  price\n",
    "even though they share the same baseline model and information\n",
    "\n",
    "- firm 1 thinks that total output will be higher and price lower than does firm 2  \n",
    "- this leads firm 1 to produce less than firm 2  \n",
    "\n",
    "\n",
    "These beliefs justify (or **rationalize**) the Markov perfect equilibrium  robust decision rules.\n",
    "\n",
    "This means that the robust rules are the unique **optimal** rules (or best responses) to the indicated worst-case transition dynamics.\n",
    "\n",
    "([[Hansen and Sargent, 2008](https://python-advanced.quantecon.org/zreferences.html#id158)] discuss how this property of robust decision rules is connected to the concept of *admissibility* in Bayesian statistical decision theory)"
   ]
  }
 ],
 "metadata": {
  "date": 1733445467.2645395,
  "filename": "rob_markov_perf.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Robust Markov Perfect Equilibrium"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}