chang_credible.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='chang-credible'></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": [
    "# Credible Government Policies in Chang Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Credible Government Policies in Chang Model](#Credible-Government-Policies-in-Chang-Model)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [The Setting](#The-Setting)  \n",
    "  - [Calculating the Set of Sustainable Promise-Value Pairs](#Calculating-the-Set-of-Sustainable-Promise-Value-Pairs)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Co-author: Sebastian Graves**\n",
    "\n",
    "In addition to what’s in Anaconda, this lecture will need the following libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": true
   },
   "outputs": [],
   "source": [
    "!pip install polytope"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "Some of the material in this lecture and [competitive equilibria in the Chang model](https://lectures.quantecon.org/chang_ramsey.html)\n",
    "can be viewed as more sophisticated and complete treatments of the topics discussed in\n",
    "[Ramsey plans, time inconsistency, sustainable plans](https://lectures.quantecon.org/calvo.html).\n",
    "\n",
    "This lecture assumes almost  the same economic environment analyzed in\n",
    "[competitive equilibria in the Chang model](https://lectures.quantecon.org/chang_ramsey.html).\n",
    "\n",
    "The only change  –  and it is a substantial one – is the timing protocol for making government decisions.\n",
    "\n",
    "In  [competitive equilibria in the Chang model](https://lectures.quantecon.org/chang_ramsey.html), a *Ramsey planner*\n",
    "chose a comprehensive government policy once-and-for-all at time $ 0 $.\n",
    "\n",
    "Now in this lecture, there is no time $ 0 $ Ramsey planner.\n",
    "\n",
    "Instead there is a sequence of government decision-makers, one for each $ t $.\n",
    "\n",
    "The time $ t $ government decision-maker choose time $ t $ government\n",
    "actions after forecasting what future governments will do.\n",
    "\n",
    "We use the notion of a *sustainable plan* proposed in [[CK90]](https://lectures.quantecon.org/zreferences.html#chari1990sustainable),\n",
    "also referred to as a *credible public policy* in [[Sto89]](https://lectures.quantecon.org/zreferences.html#stokey1989reputation).\n",
    "\n",
    "Technically, this lecture starts where lecture\n",
    "[competitive equilibria in the Chang model](https://lectures.quantecon.org/chang_ramsey.html) on Ramsey plans\n",
    "within the Chang [[Cha98]](https://lectures.quantecon.org/zreferences.html#chang1998credible) model stopped.\n",
    "\n",
    "That lecture presents recursive representations of  *competitive equilibria* and a *Ramsey plan* for a\n",
    "version of a model of Calvo [[Cal78]](https://lectures.quantecon.org/zreferences.html#calvo1978) that Chang used to analyze and illustrate these concepts.\n",
    "\n",
    "We used two operators to characterize competitive equilibria and a Ramsey plan,\n",
    "respectively.\n",
    "\n",
    "In this lecture, we define a *credible public policy* or *sustainable plan*.\n",
    "\n",
    "Starting from a large enough initial set $ Z_0 $, we use iterations on\n",
    "Chang’s set-to-set  operator $ \\tilde D(Z) $ to\n",
    "compute a set of values associated with sustainable plans.\n",
    "\n",
    "Chang’s operator $ \\tilde D(Z) $ is closely connected with the operator\n",
    "$ D(Z) $ introduced in lecture [competitive equilibria in the Chang model](https://lectures.quantecon.org/chang_ramsey.html).\n",
    "\n",
    "- $ \\tilde D(Z) $ incorporates all of the restrictions imposed in\n",
    "  constructing the operator $ D(Z) $, but $ \\ldots $.  \n",
    "- It adds some additional restrictions  \n",
    "  \n",
    "  - these additional restrictions incorporate the idea that a plan must be *sustainable*.  \n",
    "  - *sustainable* means that the government wants to implement it at all times after all histories.  \n",
    "  \n",
    "\n",
    "\n",
    "Let’s start with some standard imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import quantecon as qe\n",
    "import polytope\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Setting\n",
    "\n",
    "We begin by reviewing the set up deployed in  [competitive equilibria in the Chang model](https://lectures.quantecon.org/chang_ramsey.html).\n",
    "\n",
    "Chang’s  model, adopted from Calvo, is designed to focus on the intertemporal trade-offs between\n",
    "the welfare benefits of deflation and the welfare costs associated with\n",
    "the high tax collections required to retire money at a rate that\n",
    "delivers deflation.\n",
    "\n",
    "A benevolent time $ 0 $ government can promote\n",
    "utility generating increases in real balances only by imposing an\n",
    "infinite sequence of sufficiently large distorting tax collections.\n",
    "\n",
    "To promote the welfare increasing effects of high real balances, the\n",
    "government wants to induce  *gradual deflation*.\n",
    "\n",
    "We start by reviewing notation.\n",
    "\n",
    "For a sequence of scalars\n",
    "$ \\vec z \\equiv \\{z_t\\}_{t=0}^\\infty $, let\n",
    "$ \\vec z^t = (z_0,  \\ldots , z_t) $,\n",
    "$ \\vec z_t = (z_t, z_{t+1}, \\ldots ) $.\n",
    "\n",
    "An infinitely lived\n",
    "representative agent and an infinitely lived government exist at dates\n",
    "$ t = 0, 1, \\ldots $.\n",
    "\n",
    "The objects in play are\n",
    "\n",
    "- an initial quantity $ M_{-1} $ of nominal money holdings  \n",
    "- a sequence of inverse money growth rates $ \\vec h $ and an associated sequence of nominal money holdings $ \\vec M $  \n",
    "- a sequence of values of money $ \\vec q $  \n",
    "- a sequence of real money holdings $ \\vec m $  \n",
    "- a sequence of total tax collections $ \\vec x $  \n",
    "- a sequence of per capita rates of consumption $ \\vec c $  \n",
    "- a sequence of per capita incomes $ \\vec y $  \n",
    "\n",
    "\n",
    "A benevolent government chooses sequences\n",
    "$ (\\vec M, \\vec h, \\vec x) $ subject to a sequence of budget\n",
    "constraints and other constraints imposed by competitive equilibrium.\n",
    "\n",
    "Given tax collection and price of money sequences, a representative household chooses\n",
    "sequences $ (\\vec c, \\vec m) $ of consumption and real balances.\n",
    "\n",
    "In competitive equilibrium, the price of money sequence $ \\vec q $ clears\n",
    "markets, thereby reconciling  decisions of the government and the\n",
    "representative household."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Household’s Problem\n",
    "\n",
    "A representative household faces a nonnegative value of money sequence\n",
    "$ \\vec q $ and sequences $ \\vec y, \\vec x $ of income and total\n",
    "tax collections, respectively.\n",
    "\n",
    "The household chooses nonnegative\n",
    "sequences $ \\vec c, \\vec M $ of consumption and nominal balances,\n",
    "respectively, to maximize\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang1'></a>\n",
    "$$\n",
    "\\sum_{t=0}^\\infty \\beta^t \\left[ u(c_t) + v(q_t M_t ) \\right] \\tag{1}\n",
    "$$\n",
    "\n",
    "subject to\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang2'></a>\n",
    "$$\n",
    "q_t M_t  \\leq y_t + q_t M_{t-1} - c_t - x_t \\tag{2}\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang3'></a>\n",
    "$$\n",
    "q_t M_t  \\leq \\bar m \\tag{3}\n",
    "$$\n",
    "\n",
    "Here $ q_t $ is the reciprocal of the price level at $ t $,\n",
    "also known as the *value of money*.\n",
    "\n",
    "Chang [[Cha98]](https://lectures.quantecon.org/zreferences.html#chang1998credible) assumes that\n",
    "\n",
    "- $ u: \\mathbb{R}_+ \\rightarrow \\mathbb{R} $ is twice continuously differentiable, strictly concave, and strictly increasing;  \n",
    "- $ v: \\mathbb{R}_+ \\rightarrow \\mathbb{R} $ is twice continuously differentiable and strictly concave;  \n",
    "- $ u'(c)_{c \\rightarrow 0}  = \\lim_{m \\rightarrow 0} v'(m) = +\\infty $;  \n",
    "- there is a finite level $ m= m^f $ such that $ v'(m^f) =0 $  \n",
    "\n",
    "\n",
    "Real balances carried out of a period equal $ m_t = q_t M_t $.\n",
    "\n",
    "Inequality [(2)](#equation-eqn-chang2) is the household’s time $ t $ budget constraint.\n",
    "\n",
    "It tells how real balances $ q_t M_t $ carried out of period $ t $ depend\n",
    "on income, consumption, taxes, and real balances $ q_t M_{t-1} $\n",
    "carried into the period.\n",
    "\n",
    "Equation [(3)](#equation-eqn-chang3) imposes an exogenous upper bound\n",
    "$ \\bar m $ on the choice of real balances, where\n",
    "$ \\bar m \\geq m^f $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Government\n",
    "\n",
    "The government chooses a sequence of inverse money growth rates with\n",
    "time $ t $ component\n",
    "$ h_t \\equiv {M_{t-1}\\over M_t} \\in \\Pi \\equiv\n",
    "[ \\underline \\pi, \\overline \\pi] $, where\n",
    "$ 0 < \\underline \\pi < 1 < { 1 \\over \\beta } \\leq \\overline \\pi $.\n",
    "\n",
    "The government faces a sequence of budget constraints with time\n",
    "$ t $ component\n",
    "\n",
    "$$\n",
    "-x_t = q_t (M_t - M_{t-1})\n",
    "$$\n",
    "\n",
    "which, by using the definitions of $ m_t $ and $ h_t $, can also\n",
    "be expressed as\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang2a'></a>\n",
    "$$\n",
    "-x_t = m_t (1-h_t) \\tag{4}\n",
    "$$\n",
    "\n",
    "The  restrictions $ m_t \\in [0, \\bar m] $ and $ h_t \\in \\Pi $ evidently\n",
    "imply that $ x_t \\in X \\equiv [(\\underline  \\pi -1)\\bar m, (\\overline \\pi -1) \\bar m] $.\n",
    "\n",
    "We define the set $ E \\equiv [0,\\bar m] \\times \\Pi \\times X $, so that we\n",
    "require that $ (m, h, x) \\in E $.\n",
    "\n",
    "To represent the idea that taxes are distorting, Chang makes the following\n",
    "assumption about outcomes for per capita output:\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang3a'></a>\n",
    "$$\n",
    "y_t = f(x_t) \\tag{5}\n",
    "$$\n",
    "\n",
    "where $ f: \\mathbb{R}\\rightarrow \\mathbb{R} $ satisfies $ f(x)  > 0 $,\n",
    "is twice continuously differentiable, $ f''(x) < 0 $, and\n",
    "$ f(x) = f(-x) $ for all $ x \\in\n",
    "\\mathbb{R} $, so that subsidies and taxes are equally distorting.\n",
    "\n",
    "The purpose is not to model the causes of tax distortions in any detail but simply to summarize\n",
    "the *outcome* of those distortions via the function $ f(x) $.\n",
    "\n",
    "A key part of the specification is that tax distortions are increasing in the\n",
    "absolute value of tax revenues.\n",
    "\n",
    "The government chooses a competitive equilibrium that\n",
    "maximizes [(1)](#equation-eqn-chang1)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Within-period Timing Protocol\n",
    "\n",
    "For the results in this lecture, the *timing* of actions within a period is\n",
    "important because of the incentives that it activates.\n",
    "\n",
    "Chang assumed the following within-period timing of decisions:\n",
    "\n",
    "- first, the government chooses $ h_t $ and $ x_t $;  \n",
    "- then given $ \\vec q $ and its expectations about future values of\n",
    "  $ x $ and $ y $’s, the household chooses $ M_t $ and therefore\n",
    "  $ m_t $ because $ m_t = q_t M_t $;  \n",
    "- then output $ y_t = f(x_t) $ is realized;  \n",
    "- finally $ c_t = y_t $  \n",
    "\n",
    "\n",
    "This within-period timing confronts the government with\n",
    "choices framed by how the private sector wants to respond when the\n",
    "government takes time $ t $ actions that differ from what the\n",
    "private sector had expected.\n",
    "\n",
    "This timing will shape the incentives confronting the government at each\n",
    "history that are to be incorporated in the construction of the $ \\tilde D $\n",
    "operator below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Household’s Problem\n",
    "\n",
    "Given $ M_{-1} $ and $ \\{q_t\\}_{t=0}^\\infty $, the household’s problem is\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathcal{L} & = \\max_{\\vec c, \\vec M}\n",
    "\\min_{\\vec \\lambda, \\vec \\mu} \\sum_{t=0}^\\infty \\beta^t\n",
    "\\bigl\\{ u(c_t) + v(M_t q_t) +\n",
    "\\lambda_t [ y_t - c_t - x_t + q_t M_{t-1} - q_t M_t ]  \\\\\n",
    "& \\quad \\quad \\quad  + \\mu_t [\\bar m - q_t  M_t] \\bigr\\}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "First-order conditions with respect to $ c_t $ and $ M_t $, respectively, are\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "u'(c_t) & = \\lambda_t \\\\\n",
    "q_t [ u'(c_t) - v'(M_t q_t) ] & \\leq \\beta u'(c_{t+1})\n",
    "q_{t+1} , \\quad = \\ {\\rm if} \\ M_t q_t < \\bar m\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Using $ h_t = {M_{t-1}\\over M_t} $ and $ q_t = {m_t \\over M_t} $ in\n",
    "these first-order conditions and rearranging implies\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang4'></a>\n",
    "$$\n",
    "m_t [u'(c_t) - v'(m_t) ] \\leq \\beta u'(f(x_{t+1})) m_{t+1} h_{t+1},\n",
    "\\quad = \\text{ if } m_t < \\bar m \\tag{6}\n",
    "$$\n",
    "\n",
    "Define the following key variable\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang5'></a>\n",
    "$$\n",
    "\\theta_{t+1} \\equiv u'(f(x_{t+1})) m_{t+1} h_{t+1} \\tag{7}\n",
    "$$\n",
    "\n",
    "This is real money balances at time $ t+1 $ measured in units of marginal\n",
    "utility, which Chang refers to as ‘the marginal utility of real balances’.\n",
    "\n",
    "From the standpoint of the household at time $ t $, equation [(7)](#equation-eqn-chang5)\n",
    "shows that $ \\theta_{t+1} $ intermediates the influences of\n",
    "$ (\\vec x_{t+1}, \\vec m_{t+1}) $ on the household’s choice of real\n",
    "balances $ m_t $.\n",
    "\n",
    "By “intermediates” we mean that the future paths\n",
    "$ (\\vec x_{t+1}, \\vec m_{t+1}) $ influence $ m_t $ entirely through\n",
    "their effects on the scalar $ \\theta_{t+1} $.\n",
    "\n",
    "The observation that the one dimensional promised marginal utility of real\n",
    "balances $ \\theta_{t+1} $ functions in this way is an important step\n",
    "in constructing a class of competitive equilibria that have a recursive representation.\n",
    "\n",
    "A closely related observation pervaded the analysis of Stackelberg plans in\n",
    "[dynamic Stackelberg problems](https://lectures.quantecon.org/dyn_stack.html) and [the Calvo model](https://lectures.quantecon.org/calvo.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Competitive Equilibrium\n",
    "\n",
    "**Definition:**\n",
    "\n",
    "- A *government policy* is a pair of sequences $ (\\vec h,\\vec x) $ where $ h_t \\in \\Pi  \\ \\forall t \\geq 0 $.  \n",
    "- A *price system* is a non-negative value of money sequence $ \\vec q $.  \n",
    "- An *allocation* is a  triple of non-negative sequences $ (\\vec c, \\vec m, \\vec y) $.  \n",
    "\n",
    "\n",
    "It is required that time $ t $ components $ (m_t, x_t, h_t) \\in E $.\n",
    "\n",
    "**Definition:**\n",
    "\n",
    "Given $ M_{-1} $, a government policy $ (\\vec h, \\vec x) $, price system $ \\vec q $, and allocation\n",
    "$ (\\vec c, \\vec m, \\vec y) $ are said to be a *competitive equilibrium* if\n",
    "\n",
    "- $ m_t = q_t M_t $ and $ y_t = f(x_t) $.  \n",
    "- The government budget constraint is satisfied.  \n",
    "- Given $ \\vec q, \\vec x, \\vec y $, $ (\\vec c, \\vec m) $ solves the household’s problem.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A Credible Government Policy\n",
    "\n",
    "Chang works with\n",
    "\n",
    "**A credible government policy with a recursive representation**\n",
    "\n",
    "- Here there is no time $ 0 $ Ramsey planner.  \n",
    "- Instead there is a sequence of governments, one for each $ t $, that\n",
    "  choose time $ t $ government actions after forecasting what future governments will do.  \n",
    "- Let $ w=\\sum_{t=0}^\\infty \\beta^t \\left[ u(c_t) + v(q_t M_t ) \\right] $\n",
    "  be a value associated with a particular competitive equilibrium.  \n",
    "- A recursive representation of a credible government policy is a pair of\n",
    "  initial conditions $ (w_0, \\theta_0) $ and a five-tuple of functions  \n",
    "  $$\n",
    "  h(w_t, \\theta_t), m(h_t, w_t, \\theta_t), x(h_t, w_t, \\theta_t), \\chi(h_t, w_t, \\theta_t),\\Psi(h_t, w_t, \\theta_t)\n",
    "  $$\n",
    "  mapping $ w_t,\\theta_t $ and in some cases $ h_t $ into\n",
    "  $ \\hat h_t, m_t, x_t, w_{t+1} $, and $ \\theta_{t+1} $, respectively.  \n",
    "- Starting from an initial condition $ (w_0, \\theta_0) $, a credible\n",
    "  government policy can be constructed by iterating on these functions in\n",
    "  the following order that respects the within-period timing:  \n",
    "  $$\n",
    "  \\begin{aligned}\n",
    "  \\hat h_t & = h(w_t,\\theta_t) \\\\\n",
    "  m_t & = m(h_t, w_t,\\theta_t) \\\\\n",
    "  x_t & = x(h_t, w_t,\\theta_t) \\\\\n",
    "  w_{t+1} & = \\chi(h_t, w_t,\\theta_t)  \\\\\n",
    "  \\theta_{t+1}  & = \\Psi(h_t, w_t,\\theta_t)\n",
    "  \\end{aligned} \\tag{8}\n",
    "  $$\n",
    "- Here it is to be understood that $ \\hat h_t $ is the action that the\n",
    "  government policy instructs the government to take, while $ h_t $\n",
    "  possibly not equal to $ \\hat h_t $ is some other action that the\n",
    "  government is free to take at time $ t $.  \n",
    "\n",
    "\n",
    "The plan is *credible* if it is in the time $ t $ government’s interest to\n",
    "execute it.\n",
    "\n",
    "Credibility requires that the plan be such that for all possible choices of\n",
    "$ h_t $ that are consistent with competitive equilibria,\n",
    "\n",
    "$$\n",
    "\\begin{split} & u(f(x(\\hat h_t, w_t,\\theta_t))) + v(m(\\hat h_t, w_t,\\theta_t))  + \\beta \\chi(\\hat h_t, w_t,\\theta_t) \\\\\n",
    "&  \\geq\n",
    "u(f(x( h_t, w_t,\\theta_t))) + v(m(h_t, w_t,\\theta_t)) + \\beta \\chi(h_t, w_t,\\theta_t) \\end{split}\n",
    "$$\n",
    "\n",
    "so that at each instance and circumstance of choice, a government attains a\n",
    "weakly higher lifetime utility with continuation value\n",
    "$ w_{t+1}=\\Psi(h_t, w_t,\\theta_t) $ by adhering to the plan and\n",
    "confirming the associated time $ t $ action $ \\hat h_t $ that\n",
    "the public had expected earlier.\n",
    "\n",
    "Please note the subtle change in arguments of the functions used to represent\n",
    "a competitive equilibrium and a Ramsey plan, on the one hand, and a credible\n",
    "government plan, on the other hand.\n",
    "\n",
    "The extra arguments appearing in the functions used to represent a credible plan\n",
    "come from allowing the government to contemplate disappointing the private sector’s\n",
    "expectation about its time $ t $ choice $ \\hat h_t $.\n",
    "\n",
    "A credible plan induces the government to confirm the private sector’s expectation.\n",
    "\n",
    "The recursive representation of the plan uses the evolution of continuation\n",
    "values to deter the government from wanting to disappoint the private sector’s\n",
    "expectations.\n",
    "\n",
    "Technically, a Ramsey plan and a credible plan  both incorporate history dependence.\n",
    "\n",
    "For a Ramsey plan, this is encoded in the dynamics of the state variable\n",
    "$ \\theta_t $, a promised marginal utility that the Ramsey plan delivers to\n",
    "the private sector.\n",
    "\n",
    "For a credible government plan, we the two-dimensional state vector\n",
    "$ (w_t, \\theta_t) $ encodes  history dependence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sustainable Plans\n",
    "\n",
    "A government strategy $ \\sigma $ and an allocation rule\n",
    "$ \\alpha $ are said to constitute a *sustainable plan* (SP) if.\n",
    "\n",
    "1. $ \\sigma $ is admissible.  \n",
    "1. Given $ \\sigma $, $ \\alpha $ is competitive.  \n",
    "1. After any history $ \\vec h^{t-1} $, the continuation of $ \\sigma $\n",
    "  is optimal for the government; i.e., the sequence $ \\vec h_t $ induced\n",
    "  by $ \\sigma $ after $ \\vec h^{t-1} $ maximizes over $ CE_\\pi $\n",
    "  given $ \\alpha $.  \n",
    "\n",
    "\n",
    "Given any history $ \\vec h^{t-1} $, the continuation of a sustainable plan is a\n",
    "sustainable plan.\n",
    "\n",
    "Let $ \\Theta = \\{ (\\vec m, \\vec x, \\vec h) \\in CE : \\text{there is an SP whose outcome is} (\\vec m, \\vec x, \\vec h) \\} $.\n",
    "\n",
    "Sustainable outcomes are elements of $ \\Theta $.\n",
    "\n",
    "Now consider the space\n",
    "\n",
    "$$\n",
    "S = \\Bigl\\{ (w,\\theta) : \\text{there is a sustainable outcome }\n",
    "    (\\vec m, \\vec x, \\vec h) \\in \\Theta\n",
    "$$\n",
    "\n",
    "with value\n",
    "\n",
    "$$\n",
    "w = \\sum_{t=0}^\\infty \\beta^t [u(f(x_t)) + v(m_t)]  \\text{ and such that }\n",
    "     u'(f(x_0)) (m_0 + x_0) = \\theta \\Bigr\\}\n",
    "$$\n",
    "\n",
    "The space $ S $ is a compact subset of $ W \\times \\Omega $\n",
    "where $ W = [\\underline w, \\overline w] $ is the space of values\n",
    "associated with sustainable plans. Here $ \\underline w $ and\n",
    "$ \\overline w $ are finite bounds on the set of values.\n",
    "\n",
    "Because there is at least one sustainable plan, $ S $ is nonempty.\n",
    "\n",
    "Now recall the within-period timing protocol, which we can depict\n",
    "$ (h,x) \\rightarrow m=q M \\rightarrow y = c $.\n",
    "\n",
    "With this timing protocol in mind, the time $ 0 $ component of an SP has the\n",
    "following components:\n",
    "\n",
    "1. A period $ 0 $ action $ \\hat h \\in \\Pi $ that the public\n",
    "  expects the government to take, together with subsequent within-period\n",
    "  consequences $ m(\\hat h), x(\\hat h) $ when the government acts as\n",
    "  expected.  \n",
    "1. For any first-period action $ h \\neq \\hat h $ with\n",
    "  $ h \\in CE_\\pi^0 $, a pair of within-period consequences\n",
    "  $ m(h), x(h) $ when the government does not act as the public had\n",
    "  expected.  \n",
    "1. For every $ h \\in \\Pi $, a pair\n",
    "  $ (w'(h), \\theta'(h))\\in S $ to carry into next period.  \n",
    "\n",
    "\n",
    "These components must be such that it is optimal for the government to\n",
    "choose $ \\hat h $ as expected; and for every possible\n",
    "$ h \\in \\Pi $, the government budget constraint and the household’s\n",
    "Euler equation must hold with continuation $ \\theta $ being\n",
    "$ \\theta'(h) $.\n",
    "\n",
    "Given the timing protocol within the model, the representative\n",
    "household’s response to a government deviation to $ h \\neq \\hat h $\n",
    "from a prescribed $ \\hat h $ consists of a first-period action\n",
    "$ m(h) $ and associated subsequent actions, together with future\n",
    "equilibrium prices, captured by $ (w'(h), \\theta'(h)) $.\n",
    "\n",
    "At this point, Chang introduces an idea in the spirit of Abreu, Pearce, and Stacchetti [[APS90]](https://lectures.quantecon.org/zreferences.html#aps1990).\n",
    "\n",
    "Let $ Z $ be a nonempty subset of $ W \\times \\Omega $.\n",
    "\n",
    "Think of using pairs $ (w', \\theta') $ drawn from $ Z $ as candidate\n",
    "continuation value, promised marginal utility pairs.\n",
    "\n",
    "Define the following operator:\n",
    "\n",
    "\n",
    "<a id='equation-chang-operator'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\tilde D(Z) = \\Bigl\\{\n",
    "(w,\\theta): \\text{there is } \\hat h \\in CE_\\pi^0 \\text{ and for each } h \\in CE_\\pi^0 \\\\\n",
    "\\text{ a four-tuple } (m(h), x(h), w'(h), \\theta'(h)) \\in [0,\\bar m] \\times X \\times Z\n",
    "\\end{aligned} \\tag{9}\n",
    "$$\n",
    "\n",
    "such that\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang12'></a>\n",
    "$$\n",
    "w = u(f(x(\\hat h)))+ v(m(\\hat h)) + \\beta w'(\\hat h) \\tag{10}\n",
    "$$\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang13'></a>\n",
    "$$\n",
    "\\theta = u'(f(x(\\hat h))) ( m(\\hat h) + x(\\hat h)) \\tag{11}\n",
    "$$\n",
    "\n",
    "and for all $ h \\in CE_\\pi^0 $\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang14'></a>\n",
    "$$\n",
    "w \\geq u(f(x(h))) + v(m(h)) + \\beta w'(h) \\tag{12}\n",
    "$$\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang-15'></a>\n",
    "$$\n",
    "x(h) = m(h) (h-1) \\tag{13}\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang16'></a>\n",
    "$$\n",
    "m(h) (u'(f(x(h))) - v'(m(h))) \\leq \\beta \\theta'(h) \\tag{14}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\quad \\quad \\ \\text{ with equality if } m(h) < \\bar m  \\Bigr\\}\n",
    "$$\n",
    "\n",
    "This operator adds the key incentive constraint to the conditions that\n",
    "had defined the earlier $ D(Z) $ operator defined in  [competitive equilibria in the Chang model](https://lectures.quantecon.org/chang_ramsey.html).\n",
    "\n",
    "Condition [(12)](#equation-eqn-chang14) requires that the plan deter the government from wanting to\n",
    "take one-shot deviations when candidate continuation values are drawn\n",
    "from $ Z $.\n",
    "\n",
    "**Proposition:**\n",
    "\n",
    "1. If $ Z \\subset \\tilde D(Z) $, then $ \\tilde D(Z) \\subset S $ (‘self-generation’).  \n",
    "1. $ S = \\tilde D(S) $ (‘factorization’).  \n",
    "\n",
    "\n",
    "**Proposition:**.\n",
    "\n",
    "1. Monotonicity of $ \\tilde D $: $ Z \\subset Z' $ implies $ \\tilde D(Z) \\subset \\tilde D(Z') $.  \n",
    "1. $ Z $ compact implies that $ \\tilde D(Z) $ is compact.  \n",
    "\n",
    "\n",
    "Chang establishes that $ S $ is compact and that therefore there\n",
    "exists a highest value SP and a lowest value SP.\n",
    "\n",
    "Further, the preceding structure allows Chang to compute $ S $ by iterating to convergence\n",
    "on $ \\tilde D $ provided that one begins with a sufficiently large\n",
    "initial set $ Z_0 $.\n",
    "\n",
    "This structure delivers the following recursive representation of a\n",
    "sustainable outcome:\n",
    "\n",
    "1. choose an initial $ (w_0, \\theta_0) \\in S $;  \n",
    "1. generate a sustainable outcome recursively by iterating on [(8)](#equation-chang501), which we repeat here for convenience:  \n",
    "\n",
    "\n",
    "> $$\n",
    "\\begin{aligned}\n",
    "\\hat h_t & = h(w_t,\\theta_t) \\\\\n",
    "m_t & = m(h_t, w_t,\\theta_t) \\\\\n",
    "x_t & = x(h_t, w_t,\\theta_t) \\\\\n",
    "w_{t+1} & = \\chi(h_t, w_t,\\theta_t)  \\\\\n",
    "\\theta_{t+1}  & = \\Psi(h_t, w_t,\\theta_t)\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating the Set of Sustainable Promise-Value Pairs\n",
    "\n",
    "Above we defined the $ \\tilde D(Z) $ operator as [(9)](#equation-chang-operator).\n",
    "\n",
    "Chang (1998) provides a method for dealing with the final three\n",
    "constraints.\n",
    "\n",
    "These incentive constraints ensure that the government wants to choose\n",
    "$ \\hat h $ as the private sector had expected it to.\n",
    "\n",
    "Chang’s simplification starts from the idea that, when considering\n",
    "whether or not to confirm the private sector’s expectation, the\n",
    "government only needs to consider the payoff of the *best* possible\n",
    "deviation.\n",
    "\n",
    "Equally, to provide incentives to the government, we only need to\n",
    "consider the harshest possible punishment.\n",
    "\n",
    "Let $ h $ denote some possible deviation. Chang defines:\n",
    "\n",
    "$$\n",
    "P(h;Z) = \\min u(f(x)) + v(m) + \\beta w'\n",
    "$$\n",
    "\n",
    "where the minimization is subject to\n",
    "\n",
    "$$\n",
    "x = m(h-1)\n",
    "$$\n",
    "\n",
    "$$\n",
    "m(h)(u'(f(x(h))) + v'(m(h))) \\leq \\beta \\theta'(h) \\text{ (with equality if } m(h) < \\bar m) \\}\n",
    "$$\n",
    "\n",
    "$$\n",
    "(m,x,w',\\theta') \\in [0,\\bar m] \\times X \\times Z\n",
    "$$\n",
    "\n",
    "For a given deviation $ h $, this problem finds the worst possible\n",
    "sustainable value.\n",
    "\n",
    "We then define:\n",
    "\n",
    "$$\n",
    "BR(Z) = \\max P(h;Z) \\text{ subject to } h \\in CE^0_\\pi\n",
    "$$\n",
    "\n",
    "$ BR(Z) $ is the value of the government’s most tempting deviation.\n",
    "\n",
    "With this in hand, we can define a new operator $ E(Z) $ that is\n",
    "equivalent to the $ \\tilde D(Z) $ operator but simpler to\n",
    "implement:\n",
    "\n",
    "$$\n",
    "E(Z) = \\Bigl\\{ (w,\\theta): \\exists  h \\in CE^0_\\pi \\text{ and } (m(h),x(h),w'(h),\\theta'(h)) \\in [0,\\bar m] \\times X \\times Z\n",
    "$$\n",
    "\n",
    "such that\n",
    "\n",
    "$$\n",
    "w = u(f(x(h))) + v(m(h)) + \\beta w'(h)\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\theta = u'(f(x(h)))(m(h) + x(h))\n",
    "$$\n",
    "\n",
    "$$\n",
    "x(h) = m(h)(h-1)\n",
    "$$\n",
    "\n",
    "$$\n",
    "m(h)(u'(f(x(h))) - v'(m(h))) \\leq \\beta \\theta'(h) \\text{ (with equality if } m(h) < \\bar m)\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "w \\geq BR(Z) \\Bigr\\}\n",
    "$$\n",
    "\n",
    "Aside from the final incentive constraint, this is the same as the\n",
    "operator in  [competitive equilibria in the Chang model](https://lectures.quantecon.org/chang_ramsey.html).\n",
    "\n",
    "Consequently, to implement this operator we just need to add one step to\n",
    "our *outer hyperplane approximation algorithm* :\n",
    "\n",
    "1. Initialize subgradients, $ H $, and hyperplane levels,\n",
    "  $ C_0 $.  \n",
    "1. Given a set of subgradients, $ H $, and hyperplane levels,\n",
    "  $ C_t $, calculate $ BR(S_t) $.  \n",
    "1. Given $ H $, $ C_t $, and $ BR(S_t) $, for each\n",
    "  subgradient $ h_i \\in H $:  \n",
    "  - Solve a linear program (described below) for each action in the\n",
    "    action space.  \n",
    "  - Find the maximum and update the corresponding hyperplane level,\n",
    "    $ C_{i,t+1} $.  \n",
    "1. If $ |C_{t+1}-C_t| > \\epsilon $, return to 2.  \n",
    "\n",
    "\n",
    "**Step 1** simply creates a large initial set $ S_0 $.\n",
    "\n",
    "Given some set $ S_t $, **Step 2** then constructs the value\n",
    "$ BR(S_t) $.\n",
    "\n",
    "To do this, we solve the following problem for each point in the action\n",
    "space $ (m_j,h_j) $:\n",
    "\n",
    "$$\n",
    "\\min_{[w',\\theta']} u(f(x_j)) + v(m_j) + \\beta w'\n",
    "$$\n",
    "\n",
    "subject to\n",
    "\n",
    "$$\n",
    "H \\cdot (w',\\theta') \\leq C_t\n",
    "$$\n",
    "\n",
    "$$\n",
    "x_j = m_j(h_j-1)\n",
    "$$\n",
    "\n",
    "$$\n",
    "m_j(u'(f(x_j)) - v'(m_j)) \\leq \\beta \\theta'\\hspace{2mm} (= \\text{if } m_j < \\bar m)\n",
    "$$\n",
    "\n",
    "This gives us a matrix of possible values, corresponding to each point\n",
    "in the action space.\n",
    "\n",
    "To find $ BR(Z) $, we minimize over the $ m $ dimension and\n",
    "maximize over the $ h $ dimension.\n",
    "\n",
    "**Step 3** then constructs the set $ S_{t+1} = E(S_t) $. The linear\n",
    "program in Step 3 is designed to construct a set $ S_{t+1} $ that is\n",
    "as large as possible while satisfying the constraints of the\n",
    "$ E(S) $ operator.\n",
    "\n",
    "To do this, for each subgradient $ h_i $, and for each point in the\n",
    "action space $ (m_j,h_j) $, we solve the following problem:\n",
    "\n",
    "$$\n",
    "\\max_{[w',\\theta']} h_i \\cdot (w,\\theta)\n",
    "$$\n",
    "\n",
    "subject to\n",
    "\n",
    "$$\n",
    "H \\cdot (w',\\theta') \\leq C_t\n",
    "$$\n",
    "\n",
    "$$\n",
    "w = u(f(x_j)) + v(m_j) + \\beta w'\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\theta = u'(f(x_j))(m_j + x_j)\n",
    "$$\n",
    "\n",
    "$$\n",
    "x_j = m_j(h_j-1)\n",
    "$$\n",
    "\n",
    "$$\n",
    "m_j(u'(f(x_j)) - v'(m_j)) \\leq \\beta \\theta'\\hspace{2mm} (= \\text{if } m_j < \\bar m)\n",
    "$$\n",
    "\n",
    "$$\n",
    "w \\geq BR(Z)\n",
    "$$\n",
    "\n",
    "This problem maximizes the hyperplane level for a given set of actions.\n",
    "\n",
    "The second part of Step 3 then finds the maximum possible hyperplane\n",
    "level across the action space.\n",
    "\n",
    "The algorithm constructs a sequence of progressively smaller sets $ S_{t+1} \\subset S_t \\subset S_{t-1} \\cdots\n",
    "\\subset S_0 $.\n",
    "\n",
    "**Step 4** ends the algorithm when the difference between these sets is\n",
    "small enough.\n",
    "\n",
    "We have created a Python class that solves the model assuming the\n",
    "following functional forms:\n",
    "\n",
    "$$\n",
    "u(c) = log(c)\n",
    "$$\n",
    "\n",
    "$$\n",
    "v(m) = \\frac{1}{500}(m \\bar m - 0.5m^2)^{0.5}\n",
    "$$\n",
    "\n",
    "$$\n",
    "f(x) = 180 - (0.4x)^2\n",
    "$$\n",
    "\n",
    "The remaining parameters $ \\{\\beta, \\bar m, \\underline h, \\bar h\\} $\n",
    "are then variables to be specified for an instance of the Chang class.\n",
    "\n",
    "Below we use the class to solve the model and plot the resulting\n",
    "equilibrium set, once with $ \\beta = 0.3 $ and once with\n",
    "$ \\beta = 0.8 $. We also plot the (larger) competitive equilibrium\n",
    "sets, which we described in  [competitive equilibria in the Chang model](https://lectures.quantecon.org/chang_ramsey.html).\n",
    "\n",
    "(We have set the number of subgradients to 10 in order to speed up the\n",
    "code for now. We can increase accuracy by increasing the number of subgradients)\n",
    "\n",
    "The following code computes sustainable plans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false,
    "html-class": "collapse"
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Author: Sebastian Graves\n",
    "\n",
    "Provides a class called ChangModel to solve different\n",
    "parameterizations of the Chang (1998) model.\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "import quantecon as qe\n",
    "import time\n",
    "\n",
    "from scipy.spatial import ConvexHull\n",
    "from scipy.optimize import linprog, minimize, minimize_scalar\n",
    "from scipy.interpolate import UnivariateSpline\n",
    "import numpy.polynomial.chebyshev as cheb\n",
    "\n",
    "\n",
    "class ChangModel:\n",
    "    \"\"\"\n",
    "    Class to solve for the competitive and sustainable sets in the Chang (1998)\n",
    "    model, for different parameterizations.\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, β, mbar, h_min, h_max, n_h, n_m, N_g):\n",
    "        # Record parameters\n",
    "        self.β, self.mbar, self.h_min, self.h_max = β, mbar, h_min, h_max\n",
    "        self.n_h, self.n_m, self.N_g = n_h, n_m, N_g\n",
    "\n",
    "        # Create other parameters\n",
    "        self.m_min = 1e-9\n",
    "        self.m_max = self.mbar\n",
    "        self.N_a = self.n_h*self.n_m\n",
    "\n",
    "        # Utility and production functions\n",
    "        uc = lambda c: np.log(c)\n",
    "        uc_p = lambda c: 1/c\n",
    "        v = lambda m: 1/500 * (mbar * m - 0.5 * m**2)**0.5\n",
    "        v_p = lambda m: 0.5/500 * (mbar * m - 0.5 * m**2)**(-0.5) * (mbar - m)\n",
    "        u = lambda h, m: uc(f(h, m)) + v(m)\n",
    "\n",
    "        def f(h, m):\n",
    "            x = m * (h - 1)\n",
    "            f = 180 - (0.4 * x)**2\n",
    "            return f\n",
    "\n",
    "        def θ(h, m):\n",
    "            x = m * (h - 1)\n",
    "            θ = uc_p(f(h, m)) * (m + x)\n",
    "            return θ\n",
    "\n",
    "        # Create set of possible action combinations, A\n",
    "        A1 = np.linspace(h_min, h_max, n_h).reshape(n_h, 1)\n",
    "        A2 = np.linspace(self.m_min, self.m_max, n_m).reshape(n_m, 1)\n",
    "        self.A = np.concatenate((np.kron(np.ones((n_m, 1)), A1),\n",
    "                                 np.kron(A2, np.ones((n_h, 1)))), axis=1)\n",
    "\n",
    "        # Pre-compute utility and output vectors\n",
    "        self.euler_vec = -np.multiply(self.A[:, 1], \\\n",
    "            uc_p(f(self.A[:, 0], self.A[:, 1])) - v_p(self.A[:, 1]))\n",
    "        self.u_vec = u(self.A[:, 0], self.A[:, 1])\n",
    "        self.Θ_vec = θ(self.A[:, 0], self.A[:, 1])\n",
    "        self.f_vec = f(self.A[:, 0], self.A[:, 1])\n",
    "        self.bell_vec = np.multiply(uc_p(f(self.A[:, 0],\n",
    "                                   self.A[:, 1])),\n",
    "                                   np.multiply(self.A[:, 1],\n",
    "                                   (self.A[:, 0] - 1))) \\\n",
    "                        + np.multiply(self.A[:, 1],\n",
    "                                      v_p(self.A[:, 1]))\n",
    "\n",
    "        # Find extrema of (w, θ) space for initial guess of equilibrium sets\n",
    "        p_vec = np.zeros(self.N_a)\n",
    "        w_vec = np.zeros(self.N_a)\n",
    "        for i in range(self.N_a):\n",
    "            p_vec[i] = self.Θ_vec[i]\n",
    "            w_vec[i] = self.u_vec[i]/(1 - β)\n",
    "\n",
    "        w_space = np.array([min(w_vec[~np.isinf(w_vec)]),\n",
    "                            max(w_vec[~np.isinf(w_vec)])])\n",
    "        p_space = np.array([0, max(p_vec[~np.isinf(w_vec)])])\n",
    "        self.p_space = p_space\n",
    "\n",
    "        # Set up hyperplane levels and gradients for iterations\n",
    "        def SG_H_V(N, w_space, p_space):\n",
    "            \"\"\"\n",
    "            This function  initializes the subgradients, hyperplane levels,\n",
    "            and extreme points of the value set by choosing an appropriate\n",
    "            origin and radius. It is based on a similar function in QuantEcon's\n",
    "            Games.jl\n",
    "            \"\"\"\n",
    "\n",
    "            # First, create a unit circle. Want points placed on [0, 2π]\n",
    "            inc = 2 * np.pi / N\n",
    "            degrees = np.arange(0, 2 * np.pi, inc)\n",
    "\n",
    "            # Points on circle\n",
    "            H = np.zeros((N, 2))\n",
    "            for i in range(N):\n",
    "                x = degrees[i]\n",
    "                H[i, 0] = np.cos(x)\n",
    "                H[i, 1] = np.sin(x)\n",
    "\n",
    "            # Then calculate origin and radius\n",
    "            o = np.array([np.mean(w_space), np.mean(p_space)])\n",
    "            r1 = max((max(w_space) - o[0])**2, (o[0] - min(w_space))**2)\n",
    "            r2 = max((max(p_space) - o[1])**2, (o[1] - min(p_space))**2)\n",
    "            r = np.sqrt(r1 + r2)\n",
    "\n",
    "            # Now calculate vertices\n",
    "            Z = np.zeros((2, N))\n",
    "            for i in range(N):\n",
    "                Z[0, i] = o[0] + r*H.T[0, i]\n",
    "                Z[1, i] = o[1] + r*H.T[1, i]\n",
    "\n",
    "            # Corresponding hyperplane levels\n",
    "            C = np.zeros(N)\n",
    "            for i in range(N):\n",
    "                C[i] = np.dot(Z[:, i], H[i, :])\n",
    "\n",
    "            return C, H, Z\n",
    "\n",
    "        C, self.H, Z = SG_H_V(N_g, w_space, p_space)\n",
    "        C = C.reshape(N_g, 1)\n",
    "        self.c0_c, self.c0_s, self.c1_c, self.c1_s = np.copy(C), np.copy(C), \\\n",
    "            np.copy(C), np.copy(C)\n",
    "        self.z0_s, self.z0_c, self.z1_s, self.z1_c = np.copy(Z), np.copy(Z), \\\n",
    "            np.copy(Z), np.copy(Z)\n",
    "\n",
    "        self.w_bnds_s, self.w_bnds_c = (w_space[0], w_space[1]), \\\n",
    "            (w_space[0], w_space[1])\n",
    "        self.p_bnds_s, self.p_bnds_c = (p_space[0], p_space[1]), \\\n",
    "            (p_space[0], p_space[1])\n",
    "\n",
    "        # Create dictionaries to save equilibrium set for each iteration\n",
    "        self.c_dic_s, self.c_dic_c = {}, {}\n",
    "        self.c_dic_s[0], self.c_dic_c[0] = self.c0_s, self.c0_c\n",
    "\n",
    "    def solve_worst_spe(self):\n",
    "        \"\"\"\n",
    "        Method to solve for BR(Z). See p.449 of Chang (1998)\n",
    "        \"\"\"\n",
    "\n",
    "        p_vec = np.full(self.N_a, np.nan)\n",
    "        c = [1, 0]\n",
    "\n",
    "        # Pre-compute constraints\n",
    "        aineq_mbar = np.vstack((self.H, np.array([0, -self.β])))\n",
    "        bineq_mbar = np.vstack((self.c0_s, 0))\n",
    "\n",
    "        aineq = self.H\n",
    "        bineq = self.c0_s\n",
    "        aeq = [[0, -self.β]]\n",
    "\n",
    "        for j in range(self.N_a):\n",
    "            # Only try if consumption is possible\n",
    "            if self.f_vec[j] > 0:\n",
    "                # If m = mbar, use inequality constraint\n",
    "                if self.A[j, 1] == self.mbar:\n",
    "                    bineq_mbar[-1] = self.euler_vec[j]\n",
    "                    res = linprog(c, A_ub=aineq_mbar, b_ub=bineq_mbar,\n",
    "                                  bounds=(self.w_bnds_s, self.p_bnds_s))\n",
    "                else:\n",
    "                    beq = self.euler_vec[j]\n",
    "                    res = linprog(c, A_ub=aineq, b_ub=bineq, A_eq=aeq, b_eq=beq,\n",
    "                                  bounds=(self.w_bnds_s, self.p_bnds_s))\n",
    "                if res.status == 0:\n",
    "                    p_vec[j] = self.u_vec[j] + self.β * res.x[0]\n",
    "\n",
    "        # Max over h and min over other variables (see Chang (1998) p.449)\n",
    "        self.br_z = np.nanmax(np.nanmin(p_vec.reshape(self.n_m, self.n_h), 0))\n",
    "\n",
    "    def solve_subgradient(self):\n",
    "        \"\"\"\n",
    "        Method to solve for E(Z). See p.449 of Chang (1998)\n",
    "        \"\"\"\n",
    "\n",
    "        # Pre-compute constraints\n",
    "        aineq_C_mbar = np.vstack((self.H, np.array([0, -self.β])))\n",
    "        bineq_C_mbar = np.vstack((self.c0_c, 0))\n",
    "\n",
    "        aineq_C = self.H\n",
    "        bineq_C = self.c0_c\n",
    "        aeq_C = [[0, -self.β]]\n",
    "\n",
    "        aineq_S_mbar = np.vstack((np.vstack((self.H, np.array([0, -self.β]))),\n",
    "                                  np.array([-self.β, 0])))\n",
    "        bineq_S_mbar = np.vstack((self.c0_s, np.zeros((2, 1))))\n",
    "\n",
    "        aineq_S = np.vstack((self.H, np.array([-self.β, 0])))\n",
    "        bineq_S = np.vstack((self.c0_s, 0))\n",
    "        aeq_S = [[0, -self.β]]\n",
    "\n",
    "        # Update maximal hyperplane level\n",
    "        for i in range(self.N_g):\n",
    "            c_a1a2_c, t_a1a2_c = np.full(self.N_a, -np.inf), \\\n",
    "                np.zeros((self.N_a, 2))\n",
    "            c_a1a2_s, t_a1a2_s = np.full(self.N_a, -np.inf), \\\n",
    "                np.zeros((self.N_a, 2))\n",
    "\n",
    "            c = [-self.H[i, 0], -self.H[i, 1]]\n",
    "\n",
    "            for j in range(self.N_a):\n",
    "                # Only try if consumption is possible\n",
    "                if self.f_vec[j] > 0:\n",
    "\n",
    "                    # COMPETITIVE EQUILIBRIA\n",
    "                    # If m = mbar, use inequality constraint\n",
    "                    if self.A[j, 1] == self.mbar:\n",
    "                        bineq_C_mbar[-1] = self.euler_vec[j]\n",
    "                        res = linprog(c, A_ub=aineq_C_mbar, b_ub=bineq_C_mbar,\n",
    "                                      bounds=(self.w_bnds_c, self.p_bnds_c))\n",
    "                    # If m < mbar, use equality constraint\n",
    "                    else:\n",
    "                        beq_C = self.euler_vec[j]\n",
    "                        res = linprog(c, A_ub=aineq_C, b_ub=bineq_C, A_eq = aeq_C,\n",
    "                                      b_eq = beq_C, bounds=(self.w_bnds_c, \\\n",
    "                                          self.p_bnds_c))\n",
    "                    if res.status == 0:\n",
    "                        c_a1a2_c[j] = self.H[i, 0] * (self.u_vec[j] \\\n",
    "                            + self.β * res.x[0]) + self.H[i, 1] * self.Θ_vec[j]\n",
    "                        t_a1a2_c[j] = res.x\n",
    "\n",
    "                    # SUSTAINABLE EQUILIBRIA\n",
    "                    # If m = mbar, use inequality constraint\n",
    "                    if self.A[j, 1] == self.mbar:\n",
    "                        bineq_S_mbar[-2] = self.euler_vec[j]\n",
    "                        bineq_S_mbar[-1] = self.u_vec[j] - self.br_z\n",
    "                        res = linprog(c, A_ub=aineq_S_mbar, b_ub=bineq_S_mbar,\n",
    "                                      bounds=(self.w_bnds_s, self.p_bnds_s))\n",
    "                    # If m < mbar, use equality constraint\n",
    "                    else:\n",
    "                        bineq_S[-1] = self.u_vec[j] - self.br_z\n",
    "                        beq_S = self.euler_vec[j]\n",
    "                        res = linprog(c, A_ub=aineq_S, b_ub=bineq_S, A_eq = aeq_S,\n",
    "                                      b_eq = beq_S, bounds=(self.w_bnds_s, \\\n",
    "                                          self.p_bnds_s))\n",
    "                    if res.status == 0:\n",
    "                        c_a1a2_s[j] = self.H[i, 0] * (self.u_vec[j] \\\n",
    "                            + self.β*res.x[0]) + self.H[i, 1] * self.Θ_vec[j]\n",
    "                        t_a1a2_s[j] = res.x\n",
    "\n",
    "            idx_c = np.where(c_a1a2_c == max(c_a1a2_c))[0][0]\n",
    "            self.z1_c[:, i] = np.array([self.u_vec[idx_c]\n",
    "                                        + self.β * t_a1a2_c[idx_c, 0],\n",
    "                                        self.Θ_vec[idx_c]])\n",
    "\n",
    "            idx_s = np.where(c_a1a2_s == max(c_a1a2_s))[0][0]\n",
    "            self.z1_s[:, i] = np.array([self.u_vec[idx_s]\n",
    "                                        + self.β * t_a1a2_s[idx_s, 0],\n",
    "                                        self.Θ_vec[idx_s]])\n",
    "\n",
    "        for i in range(self.N_g):\n",
    "            self.c1_c[i] = np.dot(self.z1_c[:, i], self.H[i, :])\n",
    "            self.c1_s[i] = np.dot(self.z1_s[:, i], self.H[i, :])\n",
    "\n",
    "    def solve_sustainable(self, tol=1e-5, max_iter=250):\n",
    "        \"\"\"\n",
    "        Method to solve for the competitive and sustainable equilibrium sets.\n",
    "        \"\"\"\n",
    "\n",
    "        t = time.time()\n",
    "        diff = tol + 1\n",
    "        iters = 0\n",
    "\n",
    "        print('### --------------- ###')\n",
    "        print('Solving Chang Model Using Outer Hyperplane Approximation')\n",
    "        print('### --------------- ### \\n')\n",
    "\n",
    "        print('Maximum difference when updating hyperplane levels:')\n",
    "\n",
    "        while diff > tol and iters < max_iter:\n",
    "            iters = iters + 1\n",
    "            self.solve_worst_spe()\n",
    "            self.solve_subgradient()\n",
    "            diff = max(np.maximum(abs(self.c0_c - self.c1_c),\n",
    "                       abs(self.c0_s - self.c1_s)))\n",
    "            print(diff)\n",
    "\n",
    "            # Update hyperplane levels\n",
    "            self.c0_c, self.c0_s = np.copy(self.c1_c), np.copy(self.c1_s)\n",
    "\n",
    "            # Update bounds for w and θ\n",
    "            wmin_c, wmax_c = np.min(self.z1_c, axis=1)[0], \\\n",
    "                np.max(self.z1_c, axis=1)[0]\n",
    "            pmin_c, pmax_c = np.min(self.z1_c, axis=1)[1], \\\n",
    "                np.max(self.z1_c, axis=1)[1]\n",
    "\n",
    "            wmin_s, wmax_s = np.min(self.z1_s, axis=1)[0], \\\n",
    "                np.max(self.z1_s, axis=1)[0]\n",
    "            pmin_S, pmax_S = np.min(self.z1_s, axis=1)[1], \\\n",
    "                np.max(self.z1_s, axis=1)[1]\n",
    "\n",
    "            self.w_bnds_s, self.w_bnds_c = (wmin_s, wmax_s), (wmin_c, wmax_c)\n",
    "            self.p_bnds_s, self.p_bnds_c = (pmin_S, pmax_S), (pmin_c, pmax_c)\n",
    "\n",
    "            # Save iteration\n",
    "            self.c_dic_c[iters], self.c_dic_s[iters] = np.copy(self.c1_c), \\\n",
    "                np.copy(self.c1_s)\n",
    "            self.iters = iters\n",
    "\n",
    "        elapsed = time.time() - t\n",
    "        print('Convergence achieved after {} iterations and {} \\\n",
    "            seconds'.format(iters, round(elapsed, 2)))\n",
    "\n",
    "    def solve_bellman(self, θ_min, θ_max, order, disp=False, tol=1e-7, maxiters=100):\n",
    "        \"\"\"\n",
    "        Continuous Method to solve the Bellman equation in section 25.3\n",
    "        \"\"\"\n",
    "        mbar = self.mbar\n",
    "\n",
    "        # Utility and production functions\n",
    "        uc = lambda c: np.log(c)\n",
    "        uc_p = lambda c: 1 / c\n",
    "        v = lambda m: 1 / 500 * (mbar * m - 0.5 * m**2)**0.5\n",
    "        v_p = lambda m: 0.5/500 * (mbar*m - 0.5 * m**2)**(-0.5) * (mbar - m)\n",
    "        u = lambda h, m: uc(f(h, m)) + v(m)\n",
    "\n",
    "        def f(h, m):\n",
    "            x = m * (h - 1)\n",
    "            f = 180 - (0.4 * x)**2\n",
    "            return f\n",
    "\n",
    "        def θ(h, m):\n",
    "            x = m * (h - 1)\n",
    "            θ = uc_p(f(h, m)) * (m + x)\n",
    "            return θ\n",
    "\n",
    "        # Bounds for Maximization\n",
    "        lb1 = np.array([self.h_min, 0, θ_min])\n",
    "        ub1 = np.array([self.h_max, self.mbar - 1e-5, θ_max])\n",
    "        lb2 = np.array([self.h_min, θ_min])\n",
    "        ub2 = np.array([self.h_max, θ_max])\n",
    "\n",
    "        # Initialize Value Function coefficients\n",
    "        # Calculate roots of Chebyshev polynomial\n",
    "        k = np.linspace(order, 1, order)\n",
    "        roots = np.cos((2 * k - 1) * np.pi / (2 * order))\n",
    "        # Scale to approximation space\n",
    "        s = θ_min + (roots - -1) / 2 * (θ_max - θ_min)\n",
    "        # Create a basis matrix\n",
    "        Φ = cheb.chebvander(roots, order - 1)\n",
    "        c = np.zeros(Φ.shape[0])\n",
    "\n",
    "        # Function to minimize and constraints\n",
    "        def p_fun(x):\n",
    "            scale = -1 + 2 * (x[2] - θ_min)/(θ_max - θ_min)\n",
    "            p_fun = - (u(x[0], x[1]) \\\n",
    "                + self.β * np.dot(cheb.chebvander(scale, order - 1), c))\n",
    "            return p_fun\n",
    "\n",
    "        def p_fun2(x):\n",
    "            scale = -1 + 2*(x[1] - θ_min)/(θ_max - θ_min)\n",
    "            p_fun = - (u(x[0],mbar) \\\n",
    "                + self.β * np.dot(cheb.chebvander(scale, order - 1), c))\n",
    "            return p_fun\n",
    "\n",
    "        cons1 = ({'type': 'eq',   'fun': lambda x: uc_p(f(x[0], x[1])) * x[1]\n",
    "                    * (x[0] - 1) + v_p(x[1]) * x[1] + self.β * x[2] - θ},\n",
    "                 {'type': 'eq',   'fun': lambda x: uc_p(f(x[0], x[1]))\n",
    "                    * x[0] * x[1] - θ})\n",
    "        cons2 = ({'type': 'ineq', 'fun': lambda x: uc_p(f(x[0], mbar)) * mbar\n",
    "                    * (x[0] - 1) + v_p(mbar) * mbar + self.β * x[1] - θ},\n",
    "                 {'type': 'eq',   'fun': lambda x: uc_p(f(x[0], mbar))\n",
    "                    * x[0] * mbar - θ})\n",
    "\n",
    "        bnds1 = np.concatenate([lb1.reshape(3, 1), ub1.reshape(3, 1)], axis=1)\n",
    "        bnds2 = np.concatenate([lb2.reshape(2, 1), ub2.reshape(2, 1)], axis=1)\n",
    "\n",
    "        # Bellman Iterations\n",
    "        diff = 1\n",
    "        iters = 1\n",
    "\n",
    "        while diff > tol:\n",
    "        # 1. Maximization, given value function guess\n",
    "            p_iter1 = np.zeros(order)\n",
    "            for i in range(order):\n",
    "                θ = s[i]\n",
    "                res = minimize(p_fun,\n",
    "                               lb1 + (ub1-lb1) / 2,\n",
    "                               method='SLSQP',\n",
    "                               bounds=bnds1,\n",
    "                               constraints=cons1,\n",
    "                               tol=1e-10)\n",
    "                if res.success == True:\n",
    "                    p_iter1[i] = -p_fun(res.x)\n",
    "                res = minimize(p_fun2,\n",
    "                               lb2 + (ub2-lb2) / 2,\n",
    "                               method='SLSQP',\n",
    "                               bounds=bnds2,\n",
    "                               constraints=cons2,\n",
    "                               tol=1e-10)\n",
    "                if -p_fun2(res.x) > p_iter1[i] and res.success == True:\n",
    "                    p_iter1[i] = -p_fun2(res.x)\n",
    "\n",
    "            # 2. Bellman updating of Value Function coefficients\n",
    "            c1 = np.linalg.solve(Φ, p_iter1)\n",
    "            # 3. Compute distance and update\n",
    "            diff = np.linalg.norm(c - c1)\n",
    "            if bool(disp == True):\n",
    "                print(diff)\n",
    "            c = np.copy(c1)\n",
    "            iters = iters + 1\n",
    "            if iters > maxiters:\n",
    "                print('Convergence failed after {} iterations'.format(maxiters))\n",
    "                break\n",
    "\n",
    "        self.θ_grid = s\n",
    "        self.p_iter = p_iter1\n",
    "        self.Φ = Φ\n",
    "        self.c = c\n",
    "        print('Convergence achieved after {} iterations'.format(iters))\n",
    "\n",
    "        # Check residuals\n",
    "        θ_grid_fine = np.linspace(θ_min, θ_max, 100)\n",
    "        resid_grid = np.zeros(100)\n",
    "        p_grid = np.zeros(100)\n",
    "        θ_prime_grid = np.zeros(100)\n",
    "        m_grid = np.zeros(100)\n",
    "        h_grid = np.zeros(100)\n",
    "        for i in range(100):\n",
    "            θ = θ_grid_fine[i]\n",
    "            res = minimize(p_fun,\n",
    "                           lb1 + (ub1-lb1) / 2,\n",
    "                           method='SLSQP',\n",
    "                           bounds=bnds1,\n",
    "                           constraints=cons1,\n",
    "                           tol=1e-10)\n",
    "            if res.success == True:\n",
    "                p = -p_fun(res.x)\n",
    "                p_grid[i] = p\n",
    "                θ_prime_grid[i] = res.x[2]\n",
    "                h_grid[i] = res.x[0]\n",
    "                m_grid[i] = res.x[1]\n",
    "            res = minimize(p_fun2,\n",
    "                           lb2 + (ub2-lb2)/2,\n",
    "                           method='SLSQP',\n",
    "                           bounds=bnds2,\n",
    "                           constraints=cons2,\n",
    "                           tol=1e-10)\n",
    "            if -p_fun2(res.x) > p and res.success == True:\n",
    "                p = -p_fun2(res.x)\n",
    "                p_grid[i] = p\n",
    "                θ_prime_grid[i] = res.x[1]\n",
    "                h_grid[i] = res.x[0]\n",
    "                m_grid[i] = self.mbar\n",
    "            scale = -1 + 2 * (θ - θ_min)/(θ_max - θ_min)\n",
    "            resid_grid[i] = np.dot(cheb.chebvander(scale, order-1), c) - p\n",
    "\n",
    "        self.resid_grid = resid_grid\n",
    "        self.θ_grid_fine = θ_grid_fine\n",
    "        self.θ_prime_grid = θ_prime_grid\n",
    "        self.m_grid = m_grid\n",
    "        self.h_grid = h_grid\n",
    "        self.p_grid = p_grid\n",
    "        self.x_grid = m_grid * (h_grid - 1)\n",
    "\n",
    "        # Simulate\n",
    "        θ_series = np.zeros(31)\n",
    "        m_series = np.zeros(30)\n",
    "        h_series = np.zeros(30)\n",
    "\n",
    "        # Find initial θ\n",
    "        def ValFun(x):\n",
    "            scale = -1 + 2*(x - θ_min)/(θ_max - θ_min)\n",
    "            p_fun = np.dot(cheb.chebvander(scale, order - 1), c)\n",
    "            return -p_fun\n",
    "\n",
    "        res = minimize(ValFun,\n",
    "                      (θ_min + θ_max)/2,\n",
    "                      bounds=[(θ_min, θ_max)])\n",
    "        θ_series[0] = res.x\n",
    "\n",
    "        # Simulate\n",
    "        for i in range(30):\n",
    "            θ = θ_series[i]\n",
    "            res = minimize(p_fun,\n",
    "                           lb1 + (ub1-lb1)/2,\n",
    "                           method='SLSQP',\n",
    "                           bounds=bnds1,\n",
    "                           constraints=cons1,\n",
    "                           tol=1e-10)\n",
    "            if res.success == True:\n",
    "                p = -p_fun(res.x)\n",
    "                h_series[i] = res.x[0]\n",
    "                m_series[i] = res.x[1]\n",
    "                θ_series[i+1] = res.x[2]\n",
    "            res2 = minimize(p_fun2,\n",
    "                            lb2 + (ub2-lb2)/2,\n",
    "                            method='SLSQP',\n",
    "                            bounds=bnds2,\n",
    "                            constraints=cons2,\n",
    "                            tol=1e-10)\n",
    "            if -p_fun2(res2.x) > p and res2.success == True:\n",
    "                h_series[i] = res2.x[0]\n",
    "                m_series[i] = self.mbar\n",
    "                θ_series[i+1] = res2.x[1]\n",
    "\n",
    "        self.θ_series = θ_series\n",
    "        self.m_series = m_series\n",
    "        self.h_series = h_series\n",
    "        self.x_series = m_series * (h_series - 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparison of Sets\n",
    "\n",
    "The set of $ (w, \\theta) $ associated with  sustainable plans is  smaller than the set of $ (w, \\theta) $\n",
    "pairs associated with competitive equilibria, since the additional\n",
    "constraints associated with sustainability must also be satisfied.\n",
    "\n",
    "Let’s compute two examples, one with a low $ \\beta $, another with a higher $ \\beta $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ch1 = ChangModel(β=0.3, mbar=30, h_min=0.9, h_max=2, n_h=8, n_m=35, N_g=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ch1.solve_sustainable()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following plot shows both the set of $ w,\\theta $ pairs associated with competitive equilibria (in red)\n",
    "and the smaller set of $ w,\\theta $ pairs associated with  sustainable plans (in blue)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def plot_equilibria(ChangModel):\n",
    "    \"\"\"\n",
    "    Method to plot both equilibrium sets\n",
    "    \"\"\"\n",
    "    fig, ax = plt.subplots(figsize=(7, 5))\n",
    "\n",
    "    ax.set_xlabel('w', fontsize=16)\n",
    "    ax.set_ylabel(r\"$\\theta$\", fontsize=18)\n",
    "\n",
    "    poly_S = polytope.Polytope(ChangModel.H, ChangModel.c1_s)\n",
    "    poly_C = polytope.Polytope(ChangModel.H, ChangModel.c1_c)\n",
    "    ext_C = polytope.extreme(poly_C)\n",
    "    ext_S = polytope.extreme(poly_S)\n",
    "\n",
    "    ax.fill(ext_C[:, 0], ext_C[:, 1], 'r', zorder=-1)\n",
    "    ax.fill(ext_S[:, 0], ext_S[:, 1], 'b', zorder=0)\n",
    "\n",
    "    # Add point showing Ramsey Plan\n",
    "    idx_Ramsey = np.where(ext_C[:, 0] == max(ext_C[:, 0]))[0][0]\n",
    "    R = ext_C[idx_Ramsey, :]\n",
    "    ax.scatter(R[0], R[1], 150, 'black', 'o', zorder=1)\n",
    "    w_min = min(ext_C[:, 0])\n",
    "\n",
    "    # Label Ramsey Plan slightly to the right of the point\n",
    "    ax.annotate(\"R\", xy=(R[0], R[1]),\n",
    "                xytext=(R[0] + 0.03 * (R[0] - w_min),\n",
    "                R[1]), fontsize=18)\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "plot_equilibria(ch1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evidently, the Ramsey plan, denoted by the $ R $, is not sustainable.\n",
    "\n",
    "Let’s raise the discount factor and recompute the sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ch2 = ChangModel(β=0.8, mbar=30, h_min=0.9, h_max=1/0.8,\n",
    "    n_h=8, n_m=35, N_g=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ch2.solve_sustainable()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s plot both sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "plot_equilibria(ch2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evidently, the Ramsey plan is now sustainable."
   ]
  }
 ],
 "metadata": {
  "filename": "chang_credible.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Credible Government Policies in Chang Model"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}