chang_ramsey.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='chang-ramsey'></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": [
    "# Competitive Equilibria of a Model of Chang"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Competitive Equilibria of a Model of Chang](#Competitive-Equilibria-of-a-Model-of-Chang)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [Setting](#Setting)  \n",
    "  - [Competitive Equilibrium](#Competitive-Equilibrium)  \n",
    "  - [Inventory of Objects in Play](#Inventory-of-Objects-in-Play)  \n",
    "  - [Analysis](#Analysis)  \n",
    "  - [Calculating all Promise-Value Pairs in CE](#Calculating-all-Promise-Value-Pairs-in-CE)  \n",
    "  - [Solving a Continuation Ramsey Planner’s Bellman Equation](#Solving-a-Continuation-Ramsey-Planner’s-Bellman-Equation)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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",
    "This lecture describes how Chang [[Cha98]](https://python.quantecon.org/zreferences.html#chang1998credible)\n",
    "analyzed **competitive equilibria** and a best competitive equilibrium called a  **Ramsey plan**.\n",
    "\n",
    "He did this by\n",
    "\n",
    "- characterizing a competitive equilibrium recursively in a way also employed\n",
    "  in the [dynamic Stackelberg problems](https://python.quantecon.org/dyn_stack.html) and [Calvo model](https://python.quantecon.org/calvo.html) lectures\n",
    "  to pose Stackelberg problems in linear economies, and then  \n",
    "- appropriately adapting an argument of Abreu, Pearce, and Stachetti\n",
    "  [[APS90]](https://python.quantecon.org/zreferences.html#aps1990) to describe key features of the  set of competitive equilibria  \n",
    "\n",
    "\n",
    "Roberto Chang [[Cha98]](https://python.quantecon.org/zreferences.html#chang1998credible)  chose a model of  Calvo [[Cal78]](https://python.quantecon.org/zreferences.html#calvo1978)\n",
    "as a simple structure that conveys ideas that apply more broadly.\n",
    "\n",
    "A textbook version of Chang’s model appears in chapter 25 of [[LS18]](https://python.quantecon.org/zreferences.html#ljungqvist2012).\n",
    "\n",
    "This lecture and [Credible Government Policies in Chang Model](https://python.quantecon.org/chang_credible.html)\n",
    "can be viewed as more sophisticated and complete treatments of the\n",
    "topics discussed in [Ramsey plans, time inconsistency, sustainable plans](https://python.quantecon.org/calvo.html).\n",
    "\n",
    "Both this lecture and [Credible Government Policies in Chang Model](https://python.quantecon.org/chang_credible.html)\n",
    "make extensive use of an idea to which we apply the nickname\n",
    "**dynamic programming squared**.\n",
    "\n",
    "In dynamic programming squared problems there are typically two interrelated Bellman equations\n",
    "\n",
    "- A Bellman equation for a set of agents or followers with value or value function $ v_a $.  \n",
    "- A Bellman equation for a principal or Ramsey planner or Stackelberg leader with value or\n",
    "  value function $ v_p $ in which $ v_a $ appears as an argument.  \n",
    "\n",
    "\n",
    "We encountered problems with this structure in\n",
    "[dynamic Stackelberg problems](https://python.quantecon.org/dyn_stack.html),\n",
    "[optimal taxation with state-contingent debt](https://python.quantecon.org/opt_tax_recur.html),\n",
    "and other lectures.\n",
    "\n",
    "We’ll start with some standard imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import polytope\n",
    "import quantecon as qe\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Setting\n",
    "\n",
    "First, we introduce some 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.\n",
    "\n",
    "Chang adopts a version of a model that [[Cal78]](https://python.quantecon.org/zreferences.html#calvo1978) designed to exhibit\n",
    "time-inconsistency of a Ramsey policy in a simple and transparent\n",
    "setting.\n",
    "\n",
    "By influencing the representative household’s expectations, government actions at\n",
    "time $ t $ affect components of household utilities for periods\n",
    "$ s $ before $ t $.\n",
    "\n",
    "When setting a path for monetary\n",
    "expansion rates, the government takes into account how the\n",
    "household’s anticipations of the government’s future actions affect the household’s current decisions.\n",
    "\n",
    "The ultimate source of time inconsistency is that a\n",
    "time $ 0 $ Ramsey planner  takes these\n",
    "effects into account in designing a plan of government actions for\n",
    "$ t \\geq 0 $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting"
   ]
  },
  {
   "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-chang-ramsey1'></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-chang-ramsey2'></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-chang-ramsey3'></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",
    "which we can also call the *value of money*.\n",
    "\n",
    "Chang [[Cha98]](https://python.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",
    "The household carries real balances out of a period equal to $ m_t = q_t M_t $.\n",
    "\n",
    "Inequality [(2)](#equation-eqn-chang-ramsey2) 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-chang-ramsey3) imposes an exogenous upper bound\n",
    "$ \\bar m $ on the household’s 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-chang-ramsey2a'></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,\n",
    "(\\overline \\pi -1) \\bar m] $.\n",
    "\n",
    "We define the set $ E \\equiv [0,\\bar m] \\times \\Pi \\times X $,\n",
    "so that we 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-chang-ramsey3a'></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",
    "Calvo’s and Chang’s  purpose is not to model the causes of tax distortions in\n",
    "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",
    "**Ramsey plan:**\n",
    "A Ramsey plan is a competitive equilibrium that\n",
    "maximizes [(1)](#equation-eqn-chang-ramsey1).\n",
    "\n",
    "Within-period timing of decisions is as follows:\n",
    "\n",
    "- first, the government chooses $ h_t $ and $ x_t $;  \n",
    "- then given $ \\vec q $ and its expectations about future values of $ x $\n",
    "  and $ y $’s, the household chooses $ M_t $ and therefore $ m_t $\n",
    "  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 consideration will be important in lecture [credible government policies](https://python.quantecon.org/chang_credible.html) when\n",
    "we study *credible government policies*.\n",
    "\n",
    "The model 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 sufficiently\n",
    "large distorting tax collections.\n",
    "\n",
    "To promote the welfare increasing effects of high real balances, the\n",
    "government wants to induce  *gradual deflation*."
   ]
  },
  {
   "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",
    "The last equation expresses Karush-Kuhn-Tucker complementary slackness\n",
    "conditions (see [here](https://en.wikipedia.org/wiki/Karush%E2%80%93Kuhn%E2%80%93Tucker_conditions)).\n",
    "\n",
    "These insist that the inequality is an equality at an interior solution for $ M_t $.\n",
    "\n",
    "Using $ h_t = {M_{t-1}\\over M_t} $ and $ q_t = {m_t \\over M_t} $ in these first-order conditions and rearranging implies\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang-ramsey4'></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-chang-ramsey5'></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\n",
    "balances’.\n",
    "\n",
    "From the standpoint of the household at time $ t $, equation [(7)](#equation-eqn-chang-ramsey5)\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 balances\n",
    "$ m_t $.\n",
    "\n",
    "By “intermediates” we mean that the future paths $ (\\vec x_{t+1},\n",
    "\\vec m_{t+1}) $ influence $ m_t $ entirely through their effects on the\n",
    "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\n",
    "in lecture  [dynamic Stackelberg problems](https://python.quantecon.org/dyn_stack.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 nonnegative value of money sequence $ \\vec q $.  \n",
    "- An *allocation* is a  triple of nonnegative 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": [
    "## Inventory of Objects in Play\n",
    "\n",
    "Chang constructs the following objects\n",
    "\n",
    "1. A set $ \\Omega $ of initial marginal utilities of money $ \\theta_0 $  \n",
    "  - Let $ \\Omega $ denote the set of initial promised marginal utilities of\n",
    "    money $ \\theta_0 $ associated with competitive equilibria.  \n",
    "  - Chang exploits the fact that a competitive equilibrium consists of a first\n",
    "    period outcome $ (h_0, m_0, x_0) $ and a continuation competitive\n",
    "    equilibrium with marginal utility of money $ \\theta _1 \\in \\Omega $.  \n",
    "1. Competitive equilibria that have a recursive representation  \n",
    "  - A competitive equilibrium with a recursive representation consists of an\n",
    "    initial $ \\theta_0 $ and a four-tuple of functions $ (h, m, x, \\Psi) $\n",
    "    mapping $ \\theta $ into this period’s $ (h, m, x) $ and\n",
    "    next period’s $ \\theta $, respectively.  \n",
    "  - A competitive equilibrium can be represented recursively by iterating on  \n",
    "    $$\n",
    "    \\begin{split}\n",
    "    h_t & = h(\\theta_t) \\\\\n",
    "    m_t & = m(\\theta_t) \\\\\n",
    "    x_t & = x(\\theta_t) \\\\\n",
    "    \\theta_{t+1} & = \\Psi(\\theta_t)\n",
    "    \\end{split} \\tag{8}\n",
    "    $$\n",
    "    starting from $ \\theta_0 $  \n",
    "    The range and domain of $ \\Psi(\\cdot) $ are both $ \\Omega $  \n",
    "1. A recursive representation of a Ramsey plan  \n",
    "  - A recursive representation of a Ramsey plan is a recursive\n",
    "    competitive equilibrium $ \\theta_0, (h, m, x, \\Psi) $ that, among\n",
    "    all recursive competitive equilibria, maximizes $ \\sum_{t=0}^\\infty\n",
    "    \\beta^t \\left[ u(c_t) + v(q_t M_t ) \\right] $.  \n",
    "  - The Ramsey planner chooses $ \\theta_0, (h, m, x, \\Psi) $ from among\n",
    "    the set of recursive competitive equilibria at time $ 0 $.  \n",
    "  - Iterations on the function $ \\Psi $ determine subsequent\n",
    "    $ \\theta_t $’s that summarize the aspects of the continuation\n",
    "    competitive equilibria that influence the household’s decisions.  \n",
    "  - At time $ 0 $, the Ramsey planner commits to this implied sequence\n",
    "    $ \\{\\theta_t\\}_{t=0}^\\infty $ and therefore to an associated sequence\n",
    "    of continuation competitive equilibria.  \n",
    "1. A characterization of time-inconsistency of a Ramsey plan  \n",
    "  - Imagine that after a ‘revolution’ at time $ t \\geq 1 $, a new Ramsey\n",
    "    planner is  given the opportunity to ignore history and solve a brand new\n",
    "    Ramsey plan.  \n",
    "  - This new planner would want to reset the $ \\theta_t $ associated\n",
    "    with the original Ramsey plan to $ \\theta_0 $.  \n",
    "  - The incentive to reinitialize $ \\theta_t $ associated with this\n",
    "    revolution experiment indicates the time-inconsistency of the Ramsey plan.  \n",
    "  - By resetting $ \\theta $ to $ \\theta_0 $, the new planner avoids\n",
    "    the costs at time $ t $ that the original Ramsey planner must pay to\n",
    "    reap the beneficial effects that the original Ramsey plan for\n",
    "    $ s \\geq t $ had achieved via its influence on the household’s\n",
    "    decisions for $ s = 0, \\ldots, t-1 $.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis\n",
    "\n",
    "A competitive equilibrium is a triple of sequences\n",
    "$ (\\vec m, \\vec x, \\vec h)  \\in E^\\infty $ that satisfies\n",
    "[(2)](#equation-eqn-chang-ramsey2), [(3)](#equation-eqn-chang-ramsey3),  and [(6)](#equation-eqn-chang-ramsey4).\n",
    "\n",
    "Chang works with  a set of competitive equilibria defined as follows.\n",
    "\n",
    "**Definition:**  $ CE = \\bigl\\{ (\\vec m, \\vec x, \\vec h) \\in E^\\infty $ such that\n",
    "[(2)](#equation-eqn-chang-ramsey2), [(3)](#equation-eqn-chang-ramsey3),  and [(6)](#equation-eqn-chang-ramsey4)\n",
    "are satisfied $ \\bigr\\} $.\n",
    "\n",
    "$ CE $ is not empty because there exists a competitive equilibrium\n",
    "with $ h_t =1 $ for all $ t \\geq 1 $, namely, an equilibrium with a\n",
    "constant money supply and constant price level.\n",
    "\n",
    "Chang establishes that $ CE $ is also compact.\n",
    "\n",
    "Chang makes the following key observation that combines ideas of Abreu, Pearce,\n",
    "and Stacchetti [[APS90]](https://python.quantecon.org/zreferences.html#aps1990) with insights of Kydland and Prescott [[KP80]](https://python.quantecon.org/zreferences.html#kydland1980dynamic).\n",
    "\n",
    "**Proposition:**\n",
    "The continuation of a competitive equilibrium is a competitive equilibrium.\n",
    "\n",
    "That is, $ (\\vec m, \\vec x, \\vec h) \\in CE $ implies that $ (\\vec m_t,\n",
    "\\vec x_t, \\vec h_t) \\in CE \\  \\forall \\ t \\geq 1 $.\n",
    "\n",
    "(Lecture [dynamic Stackelberg problems](https://python.quantecon.org/dyn_stack.html) also used a version of this insight)\n",
    "\n",
    "We can now state that a **Ramsey problem** is to\n",
    "\n",
    "$$\n",
    "\\max_{(\\vec m, \\vec x, \\vec h) \\in E^\\infty} \\sum_{t=0}^\\infty \\beta^t \\left[ u(c_t) + v(m_t) \\right]\n",
    "$$\n",
    "\n",
    "subject to restrictions [(2)](#equation-eqn-chang-ramsey2), [(3)](#equation-eqn-chang-ramsey3), and [(6)](#equation-eqn-chang-ramsey4).\n",
    "\n",
    "Evidently, associated with any competitive equilibrium $ (m_0, x_0) $ is an\n",
    "implied value of $ \\theta_0 = u'(f(x_0))(m_0 + x_0) $.\n",
    "\n",
    "To bring out a recursive structure inherent in the Ramsey problem, Chang defines the set\n",
    "\n",
    "$$\n",
    "\\Omega = \\left\\{ \\theta \\in \\mathbb{R} \\\n",
    "\\text{ such that } \\ \\theta = u'(f(x_0)) (m_0 + x_0) \\ \\text{ for some } \\\n",
    "(\\vec m, \\vec x, \\vec h) \\in CE \\right\\}\n",
    "$$\n",
    "\n",
    "Equation [(6)](#equation-eqn-chang-ramsey4) inherits from the household’s Euler equation for\n",
    "money holdings the property that the value of $ m_0 $ consistent with\n",
    "the representative household’s choices depends on $ (\\vec h_1, \\vec m_1) $.\n",
    "\n",
    "This dependence is captured in the definition above by making $ \\Omega $\n",
    "be the set of first period values of $ \\theta_0 $ satisfying\n",
    "$ \\theta_0 = u'(f(x_0)) (m_0 + x_0) $ for first period component\n",
    "$ (m_0,h_0) $ of competitive equilibrium sequences\n",
    "$ (\\vec m, \\vec x, \\vec h) $.\n",
    "\n",
    "Chang establishes that $ \\Omega $ is a nonempty and compact subset of $ \\mathbb{R}_+ $.\n",
    "\n",
    "Next Chang advances:\n",
    "\n",
    "**Definition:** $ \\Gamma(\\theta) = \\{ (\\vec m, \\vec x, \\vec h) \\in CE |  \\theta = u'(f(x_0))(m_0 + x_0) \\} $.\n",
    "\n",
    "Thus, $ \\Gamma(\\theta) $ is the set of competitive equilibrium sequences\n",
    "$ (\\vec m, \\vec x, \\vec h) $ whose first period components\n",
    "$ (m_0, h_0) $ deliver the prescribed value $ \\theta $ for first\n",
    "period marginal utility.\n",
    "\n",
    "If we knew the sets $ \\Omega, \\Gamma(\\theta) $, we could use the following\n",
    "two-step procedure to find at least the *value* of the Ramsey\n",
    "outcome to the representative household\n",
    "\n",
    "1. Find the indirect value function $ w(\\theta) $ defined as  \n",
    "\n",
    "\n",
    "$$\n",
    "w(\\theta) = \\max_{(\\vec m, \\vec x, \\vec h) \\in \\Gamma(\\theta)} \\sum_{t=0}^\\infty \\beta^t \\left[ u(f(x_t)) + v(m_t) \\right]\n",
    "$$\n",
    "\n",
    "1. Compute the value of the Ramsey outcome by solving $ \\max_{\\theta \\in \\Omega} w(\\theta) $.  \n",
    "\n",
    "\n",
    "Thus, Chang states the following\n",
    "\n",
    "**Proposition**:\n",
    "\n",
    "$ w(\\theta) $ satisfies the Bellman equation\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang-ramsey7'></a>\n",
    "$$\n",
    "w(\\theta) = \\max_{x,m,h,\\theta'} \\bigl\\{ u(f(x)) + v(m) + \\beta w(\\theta') \\bigr\\} \\tag{9}\n",
    "$$\n",
    "\n",
    "where maximization is subject to\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang-ramsey8'></a>\n",
    "$$\n",
    "(m,x,h)  \\in E \\  {\\rm and} \\ \\theta' \\in \\Omega \\tag{10}\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang-ramsey9'></a>\n",
    "$$\n",
    "\\theta = u'(f(x)) (m+x) \\tag{11}\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang-ramsey10'></a>\n",
    "$$\n",
    "-x = m(1-h) \\tag{12}\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang-ramsey11'></a>\n",
    "$$\n",
    "m \\cdot [ u'(f(x)) - v'(m) ]  \\leq \\beta \\theta' , \\quad = \\ {\\rm if} \\ m < \\bar m \\tag{13}\n",
    "$$\n",
    "\n",
    "Before we use this proposition to recover a recursive representation of the\n",
    "Ramsey plan, note that the proposition relies on knowing the set $ \\Omega $.\n",
    "\n",
    "To find $ \\Omega $, Chang uses the insights of Kydland and Prescott\n",
    "[[KP80]](https://python.quantecon.org/zreferences.html#kydland1980dynamic) together with a method based on the\n",
    "Abreu, Pearce, and Stacchetti [[APS90]](https://python.quantecon.org/zreferences.html#aps1990) iteration to convergence on an\n",
    "operator $ B $ that maps continuation values into values.\n",
    "\n",
    "We want an operator that maps a continuation $ \\theta $ into a current $ \\theta $.\n",
    "\n",
    "Chang lets $ Q $ be a nonempty, bounded subset of $ \\mathbb{R} $.\n",
    "\n",
    "Elements of the set $ Q $ are taken to be candidate values for continuation marginal utilities.\n",
    "\n",
    "Chang defines an operator\n",
    "\n",
    "$$\n",
    "B(Q)  = \\theta \\in \\mathbb{R} \\\n",
    "\\text{ such that there is } \\\n",
    "(m,x,h, \\theta') \\in E \\times Q\n",
    "$$\n",
    "\n",
    "such that [(11)](#equation-eqn-chang-ramsey9), [(12)](#equation-eqn-chang-ramsey10),\n",
    "and [(13)](#equation-eqn-chang-ramsey11) hold.\n",
    "\n",
    "Thus, $ B(Q) $ is the set of first period $ \\theta $’s attainable with\n",
    "$ (m,x,h) \\in E $ and some $ \\theta' \\in Q $.\n",
    "\n",
    "**Proposition**:\n",
    "\n",
    "1. $ Q \\subset B(Q) $ implies $ B(Q) \\subset \\Omega $ (‘self-generation’).  \n",
    "1. $ \\Omega = B(\\Omega) $ (‘factorization’).  \n",
    "\n",
    "\n",
    "The proposition characterizes $ \\Omega $ as the largest fixed point\n",
    "of $ B $.\n",
    "\n",
    "It is easy to establish that $ B(Q) $ is a monotone\n",
    "operator.\n",
    "\n",
    "This property allows Chang to compute $ \\Omega $ as the\n",
    "limit of iterations on $ B $ provided that iterations begin from a\n",
    "sufficiently large initial set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Some Useful Notation\n",
    "\n",
    "Let $ \\vec h^t = (h_0, h_1, \\ldots, h_t) $ denote a history of\n",
    "inverse money creation rates with time $ t $ component\n",
    "$ h_t \\in \\Pi $.\n",
    "\n",
    "A *government strategy* $ \\sigma=\\{\\sigma_t\\}_{t=0}^\\infty $ is a $ \\sigma_0 \\in \\Pi $\n",
    "and for $ t \\geq 1 $ a sequence of functions\n",
    "$ \\sigma_t: \\Pi^{t-1} \\rightarrow \\Pi $.\n",
    "\n",
    "Chang restricts the\n",
    "government’s choice of strategies to the following space:\n",
    "\n",
    "$$\n",
    "CE_\\pi = \\{ {\\vec h} \\in \\Pi^\\infty: \\text{ there is some } \\\n",
    "(\\vec m, \\vec x) \\ \\text{ such that } \\ (\\vec m, \\vec x, \\vec h) \\in CE \\}\n",
    "$$\n",
    "\n",
    "In words, $ CE_\\pi $ is the set of money growth sequences consistent\n",
    "with the existence of competitive equilibria.\n",
    "\n",
    "Chang observes that $ CE_\\pi $ is nonempty and compact.\n",
    "\n",
    "**Definition**: $ \\sigma $ is said to be *admissible* if for all $ t \\geq 1 $\n",
    "and after any history $ \\vec h^{t-1} $, the continuation\n",
    "$ \\vec h_t $ implied by $ \\sigma $ belongs to $ CE_\\pi $.\n",
    "\n",
    "Admissibility of $ \\sigma $ means that anticipated policy choices\n",
    "associated with $ \\sigma $ are consistent with the existence of\n",
    "competitive equilibria after each possible subsequent history.\n",
    "\n",
    "After any history $ \\vec h^{t-1} $, admissibility restricts the government’s\n",
    "choice in period $ t $ to the set\n",
    "\n",
    "$$\n",
    "CE_\\pi^0 = \\{ h \\in \\Pi: {\\rm there \\ is } \\ \\vec h \\in CE_\\pi \\ {\\rm with } \\ h=h_0 \\}\n",
    "$$\n",
    "\n",
    "In words, $ CE_\\pi^0 $ is the set of all first period money growth\n",
    "rates $ h=h_0 $, each of which is consistent with the existence of a\n",
    "sequence of money growth rates $ \\vec h $ starting from $ h_0 $\n",
    "in the initial period and for which a competitive equilibrium exists.\n",
    "\n",
    "**Remark:** $ CE_\\pi^0 = \\{h \\in \\Pi: \\text{ there is } \\ (m,\\theta') \\in [0, \\bar m] \\times \\Omega \\ \\text{ such that } \\\n",
    "m u'[ f((h-1)m) - v'(m)]  \\leq \\beta \\theta' \\ \\text{ with equality if } \\  m < \\bar m \\} $.\n",
    "\n",
    "**Definition:**\n",
    "An *allocation rule* is a sequence of functions\n",
    "$ \\vec \\alpha = \\{\\alpha_t\\}_{t=0}^\\infty $ such that\n",
    "$ \\alpha_t: \\Pi^t \\rightarrow [0, \\bar m] \\times X $.\n",
    "\n",
    "Thus, the time $ t $ component of $ \\alpha_t(h^t) $ is a pair of functions\n",
    "$ (m_t(h^t), x_t(h^t)) $.\n",
    "\n",
    "**Definition:** Given an admissible government strategy\n",
    "$ \\sigma $, an allocation rule $ \\alpha $ is called\n",
    "*competitive* if given any history $ \\vec h^{t-1} $ and\n",
    "$ h_t \\in CE_\\pi^0 $, the continuations of $ \\sigma $ and\n",
    "$ \\alpha $ after $ (\\vec h^{t-1},h_t) $ induce a competitive\n",
    "equilibrium sequence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Another Operator\n",
    "\n",
    "At this point it is convenient to introduce another operator that can be\n",
    "used to compute a Ramsey plan.\n",
    "\n",
    "For computing a Ramsey plan, this\n",
    "operator is wasteful because it works with a state vector that is bigger\n",
    "than necessary.\n",
    "\n",
    "We introduce this operator because it helps to prepare\n",
    "the way for Chang’s operator called $ \\tilde D(Z) $ that we shall describe in lecture [credible government policies](https://python.quantecon.org/chang_credible.html).\n",
    "\n",
    "It is also useful because a fixed point of the operator to\n",
    "be defined here provides a good guess for an initial set\n",
    "from which to initiate iterations on Chang’s set-to-set operator  $ \\tilde D(Z) $\n",
    "to be described in lecture [credible government policies](https://python.quantecon.org/chang_credible.html).\n",
    "\n",
    "Let $ S $ be the set of all pairs $ (w, \\theta) $ of competitive\n",
    "equilibrium values and associated initial marginal utilities.\n",
    "\n",
    "Let $ W $ be a bounded set of *values* in $ \\mathbb{R} $.\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 continuation\n",
    "value, $ \\theta $ pairs.\n",
    "\n",
    "Define the operator\n",
    "\n",
    "> $$\n",
    "D(Z) = \\Bigl\\{ (w,\\theta): {\\rm there \\ is } \\ h \\in CE_\\pi^0\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\text{ and a four-tuple } \\ (m(h), x(h), w'(h), \\theta'(h)) \\in [0,\\bar m]\\times X \\times Z\n",
    "$$\n",
    "\n",
    "\n",
    "such that\n",
    "\n",
    "> \n",
    "<a id='equation-eqn-chang-ramsey120'></a>\n",
    "$$\n",
    "w = u(f(x( h))) + v(m( h)) + \\beta w'( h) \\tag{14}\n",
    "$$\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang-ramsey130'></a>\n",
    "$$\n",
    "\\theta = u'(f(x( h))) ( m( h) + x( h)) \\tag{15}\n",
    "$$\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang-ramsey-150'></a>\n",
    "$$\n",
    "x(h) = m(h) (h-1) \\tag{16}\n",
    "$$\n",
    "\n",
    "\n",
    "<a id='equation-eqn-chang-ramsey160'></a>\n",
    "$$\n",
    "m(h) (u'(f(x(h))) - v'(m(h))) \\leq \\beta \\theta'(h) \\tag{17}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\quad \\quad \\ \\text{ with equality if } m(h) < \\bar m \\Bigr\\}\n",
    "$$\n",
    "\n",
    "\n",
    "It is possible to establish.\n",
    "\n",
    "**Proposition:**\n",
    "\n",
    "1. If $ Z \\subset D(Z) $, then $ D(Z) \\subset S $ (‘self-generation’).  \n",
    "1. $ S = D(S) $ (‘factorization’).  \n",
    "\n",
    "\n",
    "**Proposition:**\n",
    "\n",
    "1. Monotonicity of $ D $: $ Z \\subset Z' $ implies $ D(Z) \\subset D(Z') $.  \n",
    "1. $ Z $ compact implies that $ D(Z) $ is compact.  \n",
    "\n",
    "\n",
    "It can be shown that $ S $ is compact and that therefore there\n",
    "exists a $ (w, \\theta) $ pair within this set that attains the\n",
    "highest possible value $ w $.\n",
    "\n",
    "This $ (w, \\theta) $ pair i\n",
    "associated with a Ramsey plan.\n",
    "\n",
    "Further, we can compute $ S $ by\n",
    "iterating to convergence on $ D $ provided that one begins with a\n",
    "sufficiently large initial set $ S_0 $.\n",
    "\n",
    "As a very useful by-product, the algorithm that finds the largest fixed\n",
    "point $ S = D(S) $ also produces the Ramsey plan, its value\n",
    "$ w $, and the associated competitive equilibrium."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating all Promise-Value Pairs in CE\n",
    "\n",
    "Above we have defined the $ D(Z) $ operator as:\n",
    "\n",
    "$$\n",
    "D(Z) = \\{ (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",
    "We noted that the set $ S $ can be found by iterating to convergence\n",
    "on $ D $, provided that we start with a sufficiently large initial\n",
    "set $ S_0 $.\n",
    "\n",
    "Our implementation builds on ideas [in this\n",
    "notebook](https://nbviewer.jupyter.org/github/QuantEcon/QuantEcon.notebooks/blob/master/recursive_repeated_games.ipynb).\n",
    "\n",
    "To find $ S $ we use a numerical algorithm called the *outer\n",
    "hyperplane approximation algorithm*.\n",
    "\n",
    "It was invented by  Judd, Yeltekin, Conklin [[JYC03]](https://python.quantecon.org/zreferences.html#juddyeltekinconklin2003).\n",
    "\n",
    "This algorithm constructs the smallest convex set that contains the\n",
    "fixed point of the $ D(S) $ operator.\n",
    "\n",
    "Given that we are finding the smallest convex set that contains\n",
    "$ S $, we can represent it on a computer as the intersection of a\n",
    "finite number of half-spaces.\n",
    "\n",
    "Let $ H $ be a set of subgradients, and $ C $ be a set of\n",
    "hyperplane levels.\n",
    "\n",
    "We approximate $ S $ by:\n",
    "\n",
    "$$\n",
    "\\tilde S = \\{(w,\\theta)| H \\cdot (w,\\theta) \\leq C \\}\n",
    "$$\n",
    "\n",
    "A key feature of this algorithm is that we discretize the action space,\n",
    "i.e., we create a grid of possible values for $ m $ and $ h $\n",
    "(note that $ x $ is implied by $ m $ and $ h $). This\n",
    "discretization simplifies computation of $ \\tilde S $ by allowing us\n",
    "to find it by solving a sequence of linear programs.\n",
    "\n",
    "The *outer hyperplane approximation algorithm* proceeds as follows:\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 $, for each 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 set\n",
    "$ S_{t+1} = D(S_t) $. The linear program in Step 2 is designed to\n",
    "construct a set $ S_{t+1} $ that is as large as possible while\n",
    "satisfying the constraints of the $ D(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",
    "This problem maximizes the hyperplane level for a given set of actions.\n",
    "\n",
    "The second part of Step 2 then finds the maximum possible hyperplane\n",
    "level across the action space.\n",
    "\n",
    "The algorithm constructs a sequence of progressively smaller sets\n",
    "$ S_{t+1} \\subset S_t \\subset S_{t-1} \\cdots\n",
    "\\subset S_0 $.\n",
    "\n",
    "**Step 3** 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 $.\n",
    "\n",
    "(Here 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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "\"\"\"\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": "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)\n",
    "ch1.solve_sustainable()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def plot_competitive(ChangModel):\n",
    "    \"\"\"\n",
    "    Method that only plots competitive equilibrium set\n",
    "    \"\"\"\n",
    "    poly_C = polytope.Polytope(ChangModel.H, ChangModel.c1_c)\n",
    "    ext_C = polytope.extreme(poly_C)\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",
    "    ax.fill(ext_C[:,0], ext_C[:,1], 'r', zorder=0)\n",
    "    ChangModel.min_theta = min(ext_C[:, 1])\n",
    "    ChangModel.max_theta = max(ext_C[:, 1])\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]), 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_competitive(ch1)"
   ]
  },
  {
   "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)\n",
    "ch2.solve_sustainable()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "plot_competitive(ch2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving a Continuation Ramsey Planner’s Bellman Equation\n",
    "\n",
    "In this section we solve the Bellman equation confronting a **continuation Ramsey planner**.\n",
    "\n",
    "The construction of a Ramsey plan is decomposed into a two subproblems in [Ramsey plans, time inconsistency, sustainable plans](https://python.quantecon.org/calvo.html)\n",
    "and [dynamic Stackelberg problems](https://python.quantecon.org/dyn_stack.html).\n",
    "\n",
    "- Subproblem 1 is faced by a sequence of continuation Ramsey planners at $ t \\geq 1 $.  \n",
    "- Subproblem 2 is faced by a Ramsey planner at $ t = 0 $.  \n",
    "\n",
    "\n",
    "The problem is:\n",
    "\n",
    "$$\n",
    "J(\\theta) = \\max_{m,x,h,\\theta'} u(f(x)) + v(m) + \\beta J(\\theta')\n",
    "$$\n",
    "\n",
    "subject to:\n",
    "\n",
    "$$\n",
    "\\theta \\leq u'(f(x))x + v'(m)m + \\beta \\theta'\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\theta = u'(f(x))(m + x )\n",
    "$$\n",
    "\n",
    "$$\n",
    "x = m(h-1)\n",
    "$$\n",
    "\n",
    "$$\n",
    "(m,x,h) \\in E\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\theta' \\in \\Omega\n",
    "$$\n",
    "\n",
    "To solve this Bellman equation, we must know the\n",
    "set $ \\Omega $.\n",
    "\n",
    "We have solved the Bellman equation for the two sets of parameter values\n",
    "for which we computed the equilibrium value sets above.\n",
    "\n",
    "Hence for these parameter configurations, we know the bounds of $ \\Omega $.\n",
    "\n",
    "The two sets of parameters\n",
    "differ only in the level of $ \\beta $.\n",
    "\n",
    "From the figures earlier in this lecture,  we know that when $ \\beta = 0.3 $,\n",
    "$ \\Omega = [0.0088,0.0499] $, and when $ \\beta = 0.8 $,\n",
    "$ \\Omega = [0.0395,0.2193] $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ch1 = ChangModel(β=0.3, mbar=30, h_min=0.99, h_max=1/0.3,\n",
    "                 n_h=8, n_m=35, N_g=50)\n",
    "ch2 = ChangModel(β=0.8, mbar=30, h_min=0.1, h_max=1/0.8,\n",
    "                 n_h=20, n_m=50, N_g=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ch1.solve_bellman(θ_min=0.01, θ_max=0.0499, order=30, tol=1e-6)\n",
    "ch2.solve_bellman(θ_min=0.045, θ_max=0.15, order=30, tol=1e-6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, a quick check that our approximations of the value functions are\n",
    "good.\n",
    "\n",
    "We do this by calculating the residuals between iterates on the value function on a fine grid:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "max(abs(ch1.resid_grid)), max(abs(ch2.resid_grid))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The value functions plotted below trace out the right edges of the sets\n",
    "of equilibrium values plotted above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n",
    "\n",
    "for ax, model in zip(axes, (ch1, ch2)):\n",
    "    ax.plot(model.θ_grid, model.p_iter)\n",
    "    ax.set(xlabel=r\"$\\theta$\",\n",
    "           ylabel=r\"$J(\\theta)$\",\n",
    "           title=rf\"$\\beta = {model.β}$\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next figure plots the optimal policy functions; values of\n",
    "$ \\theta',m,x,h $ for each value of the state $ \\theta $:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "for model in (ch1, ch2):\n",
    "\n",
    "    fig, axes = plt.subplots(2, 2, figsize=(12, 6), sharex=True)\n",
    "    fig.suptitle(rf\"$\\beta = {model.β}$\", fontsize=16)\n",
    "\n",
    "    plots = [model.θ_prime_grid, model.m_grid,\n",
    "             model.h_grid, model.x_grid]\n",
    "    labels = [r\"$\\theta'$\", \"$m$\", \"$h$\", \"$x$\"]\n",
    "\n",
    "    for ax, plot, label in zip(axes.flatten(), plots, labels):\n",
    "        ax.plot(model.θ_grid_fine, plot)\n",
    "        ax.set_xlabel(r\"$\\theta$\", fontsize=14)\n",
    "        ax.set_ylabel(label, fontsize=14)\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the first set of parameter values,  the value of $ \\theta' $ chosen by the Ramsey\n",
    "planner quickly hits the upper limit of $ \\Omega $.\n",
    "\n",
    "But with the second set of parameters  it converges to a value in the interior of the set.\n",
    "\n",
    "Consequently, the choice of $ \\bar \\theta $ is clearly important with the first set of parameter values.\n",
    "\n",
    "One way of seeing this is plotting $ \\theta'(\\theta) $ for each set\n",
    "of parameters.\n",
    "\n",
    "With the first set of parameter values,  this function does not intersect the\n",
    "45-degree line until $ \\bar \\theta $, whereas in the second set of parameter values, it\n",
    "intersects in the interior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(12, 4))\n",
    "\n",
    "for ax, model in zip(axes, (ch1, ch2)):\n",
    "    ax.plot(model.θ_grid_fine, model.θ_prime_grid, label=r\"$\\theta'(\\theta)$\")\n",
    "    ax.plot(model.θ_grid_fine, model.θ_grid_fine, label=r\"$\\theta$\")\n",
    "    ax.set(xlabel=r\"$\\theta$\", title=rf\"$\\beta = {model.β}$\")\n",
    "\n",
    "axes[0].legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Subproblem 2 is equivalent to the planner choosing the initial value of\n",
    "$ \\theta $ (i.e. the value which maximizes the value function).\n",
    "\n",
    "From this starting point, we can then trace out the paths for\n",
    "$ \\{\\theta_t,m_t,h_t,x_t\\}_{t=0}^\\infty $ that support this\n",
    "equilibrium.\n",
    "\n",
    "These are shown below for both sets of parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "for model in (ch1, ch2):\n",
    "\n",
    "    fig, axes = plt.subplots(2, 2, figsize=(12, 6))\n",
    "    fig.suptitle(rf\"$\\beta = {model.β}$\")\n",
    "\n",
    "    plots = [model.θ_series, model.m_series, model.h_series, model.x_series]\n",
    "    labels = [r\"$\\theta$\", \"$m$\", \"$h$\", \"$x$\"]\n",
    "\n",
    "    for ax, plot, label in zip(axes.flatten(), plots, labels):\n",
    "        ax.plot(plot)\n",
    "        ax.set(xlabel='t', ylabel=label)\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Next Steps\n",
    "\n",
    "In [Credible Government Policies in Chang Model](https://python.quantecon.org/chang_credible.html) we shall find\n",
    "a subset of competitive equilibria that are **sustainable**\n",
    "in the sense that a sequence of government administrations  that chooses\n",
    "sequentially, rather than once and for all at time $ 0 $ will choose to implement them.\n",
    "\n",
    "In the process of constructing them, we shall construct another,\n",
    "smaller set of competitive equilibria."
   ]
  }
 ],
 "metadata": {
  "date": 1584334735.1919355,
  "filename": "chang_ramsey.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Competitive Equilibria of a Model of Chang"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}