von_neumann_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1f3690d8",
   "metadata": {},
   "source": [
    "\n",
    "<a id='von-neumann-model'></a>\n",
    "\n",
    "<a id='index-0'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76e0d90d",
   "metadata": {},
   "source": [
    "# Von Neumann Growth Model (and a Generalization)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6224f14",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Von Neumann Growth Model (and a Generalization)](#Von-Neumann-Growth-Model-%28and-a-Generalization%29)  \n",
    "  - [Notation](#Notation)  \n",
    "  - [Model Ingredients and Assumptions](#Model-Ingredients-and-Assumptions)  \n",
    "  - [Dynamic Interpretation](#Dynamic-Interpretation)  \n",
    "  - [Duality](#Duality)  \n",
    "  - [Interpretation as Two-player Zero-sum Game](#Interpretation-as-Two-player-Zero-sum-Game)  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a6a2e5d",
   "metadata": {},
   "source": [
    "This lecture uses the class `Neumann` to calculate key objects of a\n",
    "linear growth model of John von Neumann [[von Neumann, 1937](https://python.quantecon.org/zreferences.html#id75)] that was generalized by\n",
    "Kemeny, Morgenstern and Thompson [[Kemeny *et al.*, 1956](https://python.quantecon.org/zreferences.html#id78)].\n",
    "\n",
    "Objects of interest are the maximal expansion rate ($ \\alpha $), the\n",
    "interest factor ($ β $), the optimal intensities ($ x $), and\n",
    "prices ($ p $).\n",
    "\n",
    "In addition to watching how the towering mind of John von Neumann\n",
    "formulated an equilibrium model of price and quantity vectors in\n",
    "balanced growth, this lecture shows how fruitfully to employ the\n",
    "following important tools:\n",
    "\n",
    "- a zero-sum two-player game  \n",
    "- linear programming  \n",
    "- the Perron-Frobenius theorem  \n",
    "\n",
    "\n",
    "We’ll begin with some imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9cfbf4eb",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import fsolve, linprog\n",
    "from textwrap import dedent\n",
    "\n",
    "np.set_printoptions(precision=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd32faa6",
   "metadata": {},
   "source": [
    "The code below provides the `Neumann` class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "227a3acc",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "class Neumann(object):\n",
    "\n",
    "    \"\"\"\n",
    "    This class describes the Generalized von Neumann growth model as it was\n",
    "    discussed in Kemeny et al. (1956, ECTA) and Gale (1960, Chapter 9.5):\n",
    "\n",
    "    Let:\n",
    "    n ... number of goods\n",
    "    m ... number of activities\n",
    "    A ... input matrix is m-by-n\n",
    "        a_{i,j} - amount of good j consumed by activity i\n",
    "    B ... output matrix is m-by-n\n",
    "        b_{i,j} - amount of good j produced by activity i\n",
    "\n",
    "    x ... intensity vector (m-vector) with non-negative entries\n",
    "        x'B - the vector of goods produced\n",
    "        x'A - the vector of goods consumed\n",
    "    p ... price vector (n-vector) with non-negative entries\n",
    "        Bp - the revenue vector for every activity\n",
    "        Ap - the cost of each activity\n",
    "\n",
    "    Both A and B have non-negative entries. Moreover, we assume that\n",
    "    (1) Assumption I (every good which is consumed is also produced):\n",
    "        for all j, b_{.,j} > 0, i.e. at least one entry is strictly positive\n",
    "    (2) Assumption II (no free lunch):\n",
    "        for all i, a_{i,.} > 0, i.e. at least one entry is strictly positive\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    A : array_like or scalar(float)\n",
    "        Part of the state transition equation.  It should be `n x n`\n",
    "    B : array_like or scalar(float)\n",
    "        Part of the state transition equation.  It should be `n x k`\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, A, B):\n",
    "\n",
    "        self.A, self.B = list(map(self.convert, (A, B)))\n",
    "        self.m, self.n = self.A.shape\n",
    "\n",
    "        # Check if (A, B) satisfy the basic assumptions\n",
    "        assert self.A.shape == self.B.shape, 'The input and output matrices \\\n",
    "              must have the same dimensions!'\n",
    "        assert (self.A >= 0).all() and (self.B >= 0).all(), 'The input and \\\n",
    "              output matrices must have only non-negative entries!'\n",
    "\n",
    "        # (1) Check whether Assumption I is satisfied:\n",
    "        if (np.sum(B, 0) <= 0).any():\n",
    "            self.AI = False\n",
    "        else:\n",
    "            self.AI = True\n",
    "\n",
    "        # (2) Check whether Assumption II is satisfied:\n",
    "        if (np.sum(A, 1) <= 0).any():\n",
    "            self.AII = False\n",
    "        else:\n",
    "            self.AII = True\n",
    "\n",
    "    def __repr__(self):\n",
    "        return self.__str__()\n",
    "\n",
    "    def __str__(self):\n",
    "\n",
    "        me = \"\"\"\n",
    "        Generalized von Neumann expanding model:\n",
    "          - number of goods          : {n}\n",
    "          - number of activities     : {m}\n",
    "\n",
    "        Assumptions:\n",
    "          - AI:  every column of B has a positive entry    : {AI}\n",
    "          - AII: every row of A has a positive entry       : {AII}\n",
    "\n",
    "        \"\"\"\n",
    "        # Irreducible                                       : {irr}\n",
    "        return dedent(me.format(n=self.n, m=self.m,\n",
    "                                AI=self.AI, AII=self.AII))\n",
    "\n",
    "    def convert(self, x):\n",
    "        \"\"\"\n",
    "        Convert array_like objects (lists of lists, floats, etc.) into\n",
    "        well-formed 2D NumPy arrays\n",
    "        \"\"\"\n",
    "        return np.atleast_2d(np.asarray(x))\n",
    "\n",
    "\n",
    "    def bounds(self):\n",
    "        \"\"\"\n",
    "        Calculate the trivial upper and lower bounds for alpha (expansion rate)\n",
    "        and beta (interest factor). See the proof of Theorem 9.8 in Gale (1960)\n",
    "        \"\"\"\n",
    "\n",
    "        n, m = self.n, self.m\n",
    "        A, B = self.A, self.B\n",
    "\n",
    "        f = lambda α: ((B - α * A) @ np.ones((n, 1))).max()\n",
    "        g = lambda β: (np.ones((1, m)) @ (B - β * A)).min()\n",
    "\n",
    "        UB = fsolve(f, 1).item()  # Upper bound for α, β\n",
    "        LB = fsolve(g, 2).item()  # Lower bound for α, β\n",
    "\n",
    "        return LB, UB\n",
    "\n",
    "\n",
    "    def zerosum(self, γ, dual=False):\n",
    "        \"\"\"\n",
    "        Given gamma, calculate the value and optimal strategies of a\n",
    "        two-player zero-sum game given by the matrix\n",
    "\n",
    "                M(gamma) = B - gamma * A\n",
    "\n",
    "        Row player maximizing, column player minimizing\n",
    "\n",
    "        Zero-sum game as an LP (primal --> α)\n",
    "\n",
    "            max (0', 1) @ (x', v)\n",
    "            subject to\n",
    "            [-M', ones(n, 1)] @ (x', v)' <= 0\n",
    "            (x', v) @ (ones(m, 1), 0) = 1\n",
    "            (x', v) >= (0', -inf)\n",
    "\n",
    "        Zero-sum game as an LP (dual --> beta)\n",
    "\n",
    "            min (0', 1) @ (p', u)\n",
    "            subject to\n",
    "            [M, -ones(m, 1)] @ (p', u)' <= 0\n",
    "            (p', u) @ (ones(n, 1), 0) = 1\n",
    "            (p', u) >= (0', -inf)\n",
    "\n",
    "        Outputs:\n",
    "        --------\n",
    "        value: scalar\n",
    "            value of the zero-sum game\n",
    "\n",
    "        strategy: vector\n",
    "            if dual = False, it is the intensity vector,\n",
    "            if dual = True, it is the price vector\n",
    "        \"\"\"\n",
    "\n",
    "        A, B, n, m = self.A, self.B, self.n, self.m\n",
    "        M = B - γ * A\n",
    "\n",
    "        if dual == False:\n",
    "            # Solve the primal LP (for details see the description)\n",
    "            # (1) Define the problem for v as a maximization (linprog minimizes)\n",
    "            c = np.hstack([np.zeros(m), -1])\n",
    "\n",
    "            # (2) Add constraints :\n",
    "            # ... non-negativity constraints\n",
    "            bounds = tuple(m * [(0, None)] + [(None, None)])\n",
    "            # ... inequality constraints\n",
    "            A_iq = np.hstack([-M.T, np.ones((n, 1))])\n",
    "            b_iq = np.zeros((n, 1))\n",
    "            # ... normalization\n",
    "            A_eq = np.hstack([np.ones(m), 0]).reshape(1, m + 1)\n",
    "            b_eq = 1\n",
    "\n",
    "            res = linprog(c, A_ub=A_iq, b_ub=b_iq, A_eq=A_eq, b_eq=b_eq,\n",
    "                          bounds=bounds)\n",
    "\n",
    "        else:\n",
    "            # Solve the dual LP (for details see the description)\n",
    "            # (1) Define the problem for v as a maximization (linprog minimizes)\n",
    "            c = np.hstack([np.zeros(n), 1])\n",
    "\n",
    "            # (2) Add constraints :\n",
    "            # ... non-negativity constraints\n",
    "            bounds = tuple(n * [(0, None)] + [(None, None)])\n",
    "            # ... inequality constraints\n",
    "            A_iq = np.hstack([M, -np.ones((m, 1))])\n",
    "            b_iq = np.zeros((m, 1))\n",
    "            # ... normalization\n",
    "            A_eq = np.hstack([np.ones(n), 0]).reshape(1, n + 1)\n",
    "            b_eq = 1\n",
    "\n",
    "            res = linprog(c, A_ub=A_iq, b_ub=b_iq, A_eq=A_eq, b_eq=b_eq,\n",
    "                          bounds=bounds)\n",
    "\n",
    "        if res.status != 0:\n",
    "            print(res.message)\n",
    "\n",
    "        # Pull out the required quantities\n",
    "        value = res.x[-1]\n",
    "        strategy = res.x[:-1]\n",
    "\n",
    "        return value, strategy\n",
    "\n",
    "\n",
    "    def expansion(self, tol=1e-8, maxit=1000):\n",
    "        \"\"\"\n",
    "        The algorithm used here is described in Hamburger-Thompson-Weil\n",
    "        (1967, ECTA). It is based on a simple bisection argument and utilizes\n",
    "        the idea that for a given γ (= α or β), the matrix \"M = B - γ * A\"\n",
    "        defines a two-player zero-sum game, where the optimal strategies are\n",
    "        the (normalized) intensity and price vector.\n",
    "\n",
    "        Outputs:\n",
    "        --------\n",
    "        alpha: scalar\n",
    "            optimal expansion rate\n",
    "        \"\"\"\n",
    "\n",
    "        LB, UB = self.bounds()\n",
    "\n",
    "        for iter in range(maxit):\n",
    "\n",
    "            γ = (LB + UB) / 2\n",
    "            ZS = self.zerosum(γ=γ)\n",
    "            V = ZS[0]     # value of the game with γ\n",
    "\n",
    "            if V >= 0:\n",
    "                LB = γ\n",
    "            else:\n",
    "                UB = γ\n",
    "\n",
    "            if abs(UB - LB) < tol:\n",
    "                γ = (UB + LB) / 2\n",
    "                x = self.zerosum(γ=γ)[1]\n",
    "                p = self.zerosum(γ=γ, dual=True)[1]\n",
    "                break\n",
    "\n",
    "        return γ, x, p\n",
    "\n",
    "    def interest(self, tol=1e-8, maxit=1000):\n",
    "        \"\"\"\n",
    "        The algorithm used here is described in Hamburger-Thompson-Weil\n",
    "        (1967, ECTA). It is based on a simple bisection argument and utilizes\n",
    "        the idea that for a given gamma (= alpha or beta),\n",
    "        the matrix \"M = B - γ * A\" defines a two-player zero-sum game,\n",
    "        where the optimal strategies are the (normalized) intensity and price\n",
    "        vector\n",
    "\n",
    "        Outputs:\n",
    "        --------\n",
    "        beta: scalar\n",
    "            optimal interest rate\n",
    "        \"\"\"\n",
    "\n",
    "        LB, UB = self.bounds()\n",
    "\n",
    "        for iter in range(maxit):\n",
    "            γ = (LB + UB) / 2\n",
    "            ZS = self.zerosum(γ=γ, dual=True)\n",
    "            V = ZS[0]\n",
    "\n",
    "            if V > 0:\n",
    "                LB = γ\n",
    "            else:\n",
    "                UB = γ\n",
    "\n",
    "            if abs(UB - LB) < tol:\n",
    "                γ = (UB + LB) / 2\n",
    "                p = self.zerosum(γ=γ, dual=True)[1]\n",
    "                x = self.zerosum(γ=γ)[1]\n",
    "                break\n",
    "\n",
    "        return γ, x, p"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4550a80",
   "metadata": {},
   "source": [
    "## Notation\n",
    "\n",
    "We use the following notation.\n",
    "\n",
    "$ \\mathbf{0} $ denotes\n",
    "a vector of zeros.\n",
    "\n",
    "We call an $ n $-vector  positive and write\n",
    "$ x\\gg \\mathbf{0} $ if $ x_i>0 $ for all $ i=1,2,\\dots,n $.\n",
    "\n",
    "We call a vector  non-negative and write $ x\\geq \\mathbf{0} $ if $ x_i\\geq 0 $ for\n",
    "all $ i=1,2,\\dots,n $.\n",
    "\n",
    "We call a vector  semi-positive and written $ x > \\mathbf{0} $ if\n",
    "$ x\\geq \\mathbf{0} $ and $ x\\neq \\mathbf{0} $.\n",
    "\n",
    "For two conformable vectors $ x $ and $ y $, $ x\\gg y $,\n",
    "$ x\\geq y $ and $ x> y $ mean $ x-y\\gg \\mathbf{0} $,\n",
    "$ x-y \\geq \\mathbf{0} $, and $ x-y > \\mathbf{0} $, respectively.\n",
    "\n",
    "We let all vectors in this lecture be column vectors; $ x^{T} $ denotes the\n",
    "transpose of $ x $ (i.e., a row vector).\n",
    "\n",
    "Let $ \\iota_n $ denote a\n",
    "column vector composed of $ n $ ones, i.e.\n",
    "$ \\iota_n = (1,1,\\dots,1)^T $.\n",
    "\n",
    "Let $ e^i $ denote a vector (of\n",
    "arbitrary size) containing zeros except for the $ i $ th position\n",
    "where it is one.\n",
    "\n",
    "We denote matrices by capital letters. For an arbitrary matrix\n",
    "$ A $, $ a_{i,j} $ represents the entry in its $ i $ th\n",
    "row and $ j $ th column.\n",
    "\n",
    "$ a_{\\cdot j} $ and $ a_{i\\cdot} $\n",
    "denote the $ j $ th column and $ i $ th row of $ A $,\n",
    "respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19e7abb9",
   "metadata": {},
   "source": [
    "## Model Ingredients and Assumptions\n",
    "\n",
    "A pair $ (A,B) $ of $ m\\times n $ non-negative matrices defines\n",
    "an economy.\n",
    "\n",
    "- $ m $ is the number of *activities* (or sectors)  \n",
    "- $ n $ is the number of *goods* (produced and/or consumed).  \n",
    "- $ A $ is called the *input matrix*; $ a_{i,j} $ denotes the\n",
    "  amount of good $ j $ consumed by activity $ i $  \n",
    "- $ B $ is called the *output matrix*; $ b_{i,j} $ represents\n",
    "  the amount of good $ j $ produced by activity $ i $  \n",
    "\n",
    "\n",
    "Two key assumptions restrict economy $ (A,B) $:\n",
    "\n",
    "- **Assumption I:** (every good that is consumed is also produced)  \n",
    "  $$\n",
    "  b_{.,j} > \\mathbf{0}\\hspace{5mm}\\forall j=1,2,\\dots,n\n",
    "  $$\n",
    "- **Assumption II:** (no free lunch)  \n",
    "  $$\n",
    "  a_{i,.} > \\mathbf{0}\\hspace{5mm}\\forall i=1,2,\\dots,m\n",
    "  $$\n",
    "\n",
    "\n",
    "A semi-positive  *intensity* $ m $-vector  $ x $ denotes  levels at which\n",
    "activities are operated.\n",
    "\n",
    "Therefore,\n",
    "\n",
    "- vector $ x^TA $ gives the total amount of *goods used in\n",
    "  production*  \n",
    "- vector $ x^TB $ gives *total outputs*  \n",
    "\n",
    "\n",
    "An economy $ (A,B) $ is said to be *productive*, if there exists a\n",
    "non-negative intensity vector $ x \\geq 0 $ such\n",
    "that $ x^T B > x^TA $.\n",
    "\n",
    "The semi-positive $ n $-vector $ p $ contains prices assigned to\n",
    "the $ n $ goods.\n",
    "\n",
    "The $ p $ vector implies *cost* and *revenue* vectors\n",
    "\n",
    "- the vector $ Ap $ tells *costs* of the vector of activities  \n",
    "- the vector $ Bp $ tells *revenues* from the vector of activities  \n",
    "\n",
    "\n",
    "Satisfaction or a property of an input-output pair $ (A,B) $ called *irreducibility*\n",
    "(or indecomposability) determines whether an economy can be decomposed\n",
    "into multiple     “sub-economies”.\n",
    "\n",
    "**Definition:** For an economy $ (A,B) $, the set of goods\n",
    "$ S\\subset \\{1,2,\\dots,n\\} $ is called an *independent subset* if\n",
    "it is possible to produce every good in $ S $ without consuming\n",
    "goods from outside $ S $. Formally, the set $ S $ is independent if\n",
    "$ \\exists T\\subset \\{1,2,\\dots,m\\} $ (a subset of activities) such\n",
    "that $ a_{i,j}=0 $ $ \\forall i\\in T $ and $ j\\in S^c $ and\n",
    "for all $ j\\in S $, $ \\exists i\\in T $ for which  $ b_{i,j}>0 $.\n",
    "The economy is **irreducible** if there are no proper independent\n",
    "subsets.\n",
    "\n",
    "We study two examples, both in  Chapter 9.6 of Gale [[Gale, 1989](https://python.quantecon.org/zreferences.html#id80)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56e4ff1c",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# (1) Irreducible (A, B) example: α_0 = β_0\n",
    "A1 = np.array([[0, 1, 0, 0],\n",
    "               [1, 0, 0, 1],\n",
    "               [0, 0, 1, 0]])\n",
    "\n",
    "B1 = np.array([[1, 0, 0, 0],\n",
    "               [0, 0, 2, 0],\n",
    "               [0, 1, 0, 1]])\n",
    "\n",
    "# (2) Reducible (A, B) example: β_0 < α_0\n",
    "A2 = np.array([[0, 1, 0, 0, 0, 0],\n",
    "               [1, 0, 1, 0, 0, 0],\n",
    "               [0, 0, 0, 1, 0, 0],\n",
    "               [0, 0, 1, 0, 0, 1],\n",
    "               [0, 0, 0, 0, 1, 0]])\n",
    "\n",
    "B2 = np.array([[1, 0, 0, 1, 0, 0],\n",
    "               [0, 1, 0, 0, 0, 0],\n",
    "               [0, 0, 1, 0, 0, 0],\n",
    "               [0, 0, 0, 0, 2, 0],\n",
    "               [0, 0, 0, 1, 0, 1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f3d711c",
   "metadata": {},
   "source": [
    "The following code sets up our first Neumann economy or `Neumann`\n",
    "instance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ecb76c2",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "n1 = Neumann(A1, B1)\n",
    "n1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "469859db",
   "metadata": {},
   "source": [
    "Here is a second instance of a Neumann economy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5635afe1",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "n2 = Neumann(A2, B2)\n",
    "n2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d95825a",
   "metadata": {},
   "source": [
    "## Dynamic Interpretation\n",
    "\n",
    "Attach a time index $ t $ to the preceding objects, regard an economy\n",
    "as a dynamic system, and study sequences\n",
    "\n",
    "$$\n",
    "\\{(A_t,B_t)\\}_{t\\geq 0}, \\hspace{1cm}\\{x_t\\}_{t\\geq 0},\\hspace{1cm} \\{p_t\\}_{t\\geq 0}\n",
    "$$\n",
    "\n",
    "An interesting special case holds the technology process constant and\n",
    "investigates the dynamics of quantities and prices only.\n",
    "\n",
    "Accordingly, in the rest of this lecture, we assume that\n",
    "$ (A_t,B_t)=(A,B) $ for all $ t\\geq 0 $.\n",
    "\n",
    "A crucial element of the dynamic interpretation involves the timing of\n",
    "production.\n",
    "\n",
    "We assume that production (consumption of inputs) takes place in period\n",
    "$ t $, while the consequent output materializes in period\n",
    "$ t+1 $, i.e., consumption of $ x_{t}^TA $ in period $ t $\n",
    "results in $ x^T_{t}B $ amounts of output in period $ t+1 $.\n",
    "\n",
    "These timing conventions imply the following feasibility condition:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "x^T_{t}B \\geq x^T_{t+1} A \\hspace{1cm}\\forall t\\geq 1\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "which asserts that no more goods can be used today than were produced\n",
    "yesterday.\n",
    "\n",
    "Accordingly, $ Ap_t $ tells the costs of production in period\n",
    "$ t $ and $ Bp_t $ tells revenues in period $ t+1 $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f1b2ed5",
   "metadata": {},
   "source": [
    "### Balanced Growth\n",
    "\n",
    "We follow John von Neumann in studying “balanced growth”.\n",
    "\n",
    "Let $ ./ $ denote an elementwise division of one vector by another and let\n",
    "$ \\alpha >0 $ be a scalar.\n",
    "\n",
    "Then *balanced growth* is a situation in which\n",
    "\n",
    "$$\n",
    "x_{t+1}./x_t = \\alpha , \\quad \\forall t \\geq 0\n",
    "$$\n",
    "\n",
    "With balanced growth, the law of motion of $ x $ is evidently $ x_{t+1}=\\alpha x_t $\n",
    "and so we can rewrite the feasibility constraint as\n",
    "\n",
    "$$\n",
    "x^T_{t}B \\geq \\alpha x^T_t A \\hspace{1cm}\\forall t\n",
    "$$\n",
    "\n",
    "In the same spirit, define $ \\beta\\in\\mathbb{R} $ as the **interest\n",
    "factor** per unit of time.\n",
    "\n",
    "We assume that it is always possible to earn a gross return equal to the\n",
    "constant interest factor $ \\beta $ by investing “outside the model”.\n",
    "\n",
    "Under this assumption about outside investment opportunities, a\n",
    "no-arbitrage condition gives rise to the following (no profit)\n",
    "restriction on the price sequence:\n",
    "\n",
    "$$\n",
    "\\beta Ap_{t} \\geq B p_{t} \\hspace{1cm}\\forall t\n",
    "$$\n",
    "\n",
    "This says that production cannot yield a return greater than that\n",
    "offered by the outside investment opportunity (here we compare values in\n",
    "period $ t+1 $).\n",
    "\n",
    "The balanced growth assumption allows us to drop time subscripts and\n",
    "conduct an analysis purely in terms of a time-invariant growth rate\n",
    "$ \\alpha $ and interest factor $ \\beta $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6648cb11",
   "metadata": {},
   "source": [
    "## Duality\n",
    "\n",
    "Two problems are connected by a remarkable dual\n",
    "relationship between  technological and valuation characteristics of\n",
    "the economy:\n",
    "\n",
    "**Definition:** The *technological expansion problem* (TEP) for the economy\n",
    "$ (A,B) $ is to find a semi-positive $ m $-vector $ x>0 $\n",
    "and a number $ \\alpha\\in\\mathbb{R} $ that satisfy\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    &\\max_{\\alpha} \\hspace{2mm} \\alpha\\\\\n",
    "    &\\text{s.t. }\\hspace{2mm}x^T B \\geq \\alpha x^T A\n",
    "    \\end{aligned}\n",
    "$$\n",
    "\n",
    "Theorem 9.3 of David Gale’s book [[Gale, 1989](https://python.quantecon.org/zreferences.html#id80)] asserts that if Assumptions I and II are\n",
    "both satisfied, then a maximum value of $ \\alpha $ exists and that it is\n",
    "positive.\n",
    "\n",
    "The maximal value is called the *technological expansion rate* and is denoted\n",
    "by $ \\alpha_0 $. The associated intensity vector $ x_0 $ is the\n",
    "*optimal intensity vector*.\n",
    "\n",
    "**Definition:** The economic expansion problem (EEP) for\n",
    "$ (A,B) $ is to find a semi-positive $ n $-vector $ p>0 $\n",
    "and a number $ \\beta\\in\\mathbb{R} $ that satisfy\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    &\\min_{\\beta} \\hspace{2mm} \\beta\\\\\n",
    "    &\\text{s.t. }\\hspace{2mm}Bp \\leq \\beta Ap\n",
    "    \\end{aligned}\n",
    "$$\n",
    "\n",
    "Assumptions I and II imply existence of a minimum value\n",
    "$ \\beta_0>0 $ called the *economic expansion rate*.\n",
    "\n",
    "The corresponding price vector $ p_0 $ is the *optimal price vector*.\n",
    "\n",
    "Because the criterion functions in the *technological expansion* problem\n",
    "and the *economical expansion problem* are both linearly homogeneous,\n",
    "the optimality of $ x_0 $ and $ p_0 $ are defined only up to a\n",
    "positive scale factor.\n",
    "\n",
    "For convenience (and to emphasize a close connection to zero-sum games),  we normalize both vectors\n",
    "$ x_0 $ and $ p_0 $ to have unit length.\n",
    "\n",
    "A standard duality argument (see Lemma 9.4. in (Gale, 1960) [[Gale, 1989](https://python.quantecon.org/zreferences.html#id80)]) implies\n",
    "that under Assumptions I and II, $ \\beta_0\\leq \\alpha_0 $.\n",
    "\n",
    "But to deduce that $ \\beta_0\\geq \\alpha_0 $,\n",
    "Assumptions I and II are not sufficient.\n",
    "\n",
    "Therefore, von Neumann [[von Neumann, 1937](https://python.quantecon.org/zreferences.html#id75)]  went on to prove the following remarkable\n",
    "“duality” result that connects TEP and EEP.\n",
    "\n",
    "**Theorem 1 (von Neumann):** If the economy $ (A,B) $ satisfies\n",
    "Assumptions I and II, then there exist\n",
    "$ \\left(\\gamma^{*}, x_0, p_0\\right) $, where\n",
    "$ \\gamma^{*}\\in[\\beta_0, \\alpha_0]\\subset\\mathbb{R} $, $ x_0>0 $\n",
    "is an $ m $-vector, $ p_0>0 $ is an $ n $-vector, and the\n",
    "following arbitrage true\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "x_0^T B &\\geq \\gamma^{* } x_0^T A \\\\\n",
    "Bp_0 &\\leq \\gamma^{* } Ap_0 \\\\\n",
    "x_0^T\\left(B-\\gamma^{* } A\\right)p_0 &= 0\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    ">**Note**\n",
    ">\n",
    ">*Proof (Sketch):* Assumption I and II imply that there exist $ (\\alpha_0,\n",
    "x_0) $ and $ (\\beta_0, p_0) $ that solve the TEP and EEP, respectively. If\n",
    "$ \\gamma^*>\\alpha_0 $, then by definition of $ \\alpha_0 $, there cannot\n",
    "exist a semi-positive $ x $ that satisfies $ x^T B \\geq \\gamma^{* }\n",
    "x^T A $.  Similarly, if $ \\gamma^*<\\beta_0 $, there is no semi-positive\n",
    "$ p $ for which $ Bp \\leq \\gamma^{* } Ap $. Let $ \\gamma^{*\n",
    "}\\in[\\beta_0, \\alpha_0] $, then $ x_0^T B \\geq \\alpha_0 x_0^T A \\geq\n",
    "\\gamma^{* } x_0^T A $.  Moreover, $ Bp_0\\leq \\beta_0 A p_0\\leq \\gamma^* A\n",
    "p_0 $. These two inequalities imply $ x_0\\left(B - \\gamma^{* } A\\right)p_0\n",
    "= 0 $.\n",
    "\n",
    "Here the constant $ \\gamma^{*} $ is both an expansion factor  and an interest\n",
    "factor (not necessarily optimal).\n",
    "\n",
    "We have already encountered and\n",
    "discussed the first two inequalities that represent feasibility and\n",
    "no-profit conditions.\n",
    "\n",
    "Moreover, the equality $ x_0^T\\left(B-\\gamma^{* } A\\right)p_0 = 0 $ concisely expresses the\n",
    "requirements that if any good grows at a rate larger than\n",
    "$ \\gamma^{*} $ (i.e., if it is *oversupplied*), then its price\n",
    "must be zero; and that if any activity provides negative profit, it must\n",
    "be unused.\n",
    "\n",
    "Therefore, the conditions stated in Theorem I ex encode all equilibrium conditions.\n",
    "\n",
    "So  Theorem I essentially states that under Assumptions I and II there\n",
    "always exists an equilibrium $ \\left(\\gamma^{*}, x_0, p_0\\right) $\n",
    "with balanced growth.\n",
    "\n",
    "Note that Theorem I is silent about uniqueness of the equilibrium. In\n",
    "fact, it does not rule out (trivial) cases with $ x_0^TBp_0 = 0 $ so\n",
    "that nothing of value is produced.\n",
    "\n",
    "To exclude such uninteresting cases,\n",
    "Kemeny, Morgenstern and Thomspson [[Kemeny *et al.*, 1956](https://python.quantecon.org/zreferences.html#id78)] add an extra requirement\n",
    "\n",
    "$$\n",
    "x^T_0 B p_0 > 0\n",
    "$$\n",
    "\n",
    "and call the associated equilibria *economic solutions*.\n",
    "\n",
    "They show that\n",
    "this extra condition does not affect the existence result, while it\n",
    "significantly reduces the number of (relevant) solutions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50df0e4e",
   "metadata": {},
   "source": [
    "## Interpretation as Two-player Zero-sum Game\n",
    "\n",
    "To compute the equilibrium $ (\\gamma^{*}, x_0, p_0) $, we follow the\n",
    "algorithm proposed by Hamburger, Thompson and Weil (1967), building on\n",
    "the key insight that an equilibrium (with balanced growth) can be\n",
    "solves a particular two-player zero-sum game.\n",
    "First, we introduce some notation.\n",
    "\n",
    "Consider the $ m\\times n $ matrix $ C $ as a payoff matrix,\n",
    "with the entries representing payoffs from the **minimizing** column\n",
    "player to the **maximizing** row player and assume that the players can\n",
    "use mixed strategies. Thus,\n",
    "\n",
    "- the  row player chooses the $ m $-vector $ x > \\mathbf{0} $ subject to $ \\iota_m^T x = 1 $  \n",
    "- the column player chooses the $ n $-vector $ p > \\mathbf{0} $ subject to $ \\iota_n^T p = 1 $.  \n",
    "\n",
    "\n",
    "**Definition:** The $ m\\times n $ matrix game $ C $ has the\n",
    "*solution* $ (x^*, p^*, V(C)) $ in mixed strategies if\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "(x^* )^T C e^j \\geq V(C)\\quad \\forall j\\in\\{1, \\dots, n\\}\\quad \\quad\n",
    "\\text{and}\\quad\\quad (e^i)^T C p^* \\leq V(C)\\quad \\forall i\\in\\{1, \\dots, m\\}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The number $ V(C) $ is called the *value* of the game.\n",
    "\n",
    "From the above definition, it is clear that the value $ V(C) $ has\n",
    "two alternative interpretations:\n",
    "\n",
    "- by playing the appropriate mixed\n",
    "  stategy, the maximizing player can assure himself at least $ V(C) $\n",
    "  (no matter what the column player chooses)  \n",
    "- by playing the appropriate\n",
    "  mixed stategy, the minimizing player can make sure that the maximizing\n",
    "  player will not get more than $ V(C) $ (irrespective of what is the\n",
    "  maximizing player’s choice)  \n",
    "\n",
    "\n",
    "A famous theorem of Nash (1951) tells us that there always\n",
    "exists a mixed strategy Nash equilibrium for any *finite* two-player\n",
    "zero-sum game.\n",
    "\n",
    "Moreover, von Neumann’s Minmax Theorem [[von Neumann, 1928](https://python.quantecon.org/zreferences.html#id76)] implies that\n",
    "\n",
    "$$\n",
    "V(C) = \\max_x \\min_p \\hspace{2mm} x^T C p = \\min_p \\max_x \\hspace{2mm} x^T C p = (x^*)^T C p^*\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "791fda25",
   "metadata": {},
   "source": [
    "### Connection with Linear Programming (LP)\n",
    "\n",
    "Nash equilibria of a finite two-player zero-sum game solve  a linear programming problem.\n",
    "\n",
    "To see this, we introduce\n",
    "the following notation\n",
    "\n",
    "- For a fixed $ x $, let $ v $ be the value of the minimization problem: $ v \\equiv \\min_p x^T C p = \\min_j x^T C e^j $  \n",
    "- For a fixed $ p $, let $ u $ be the value of the maximization problem: $ u \\equiv \\max_x x^T C p = \\max_i (e^i)^T C p $  \n",
    "\n",
    "\n",
    "Then the *max-min problem* (the game from the maximizing player’s point\n",
    "of view) can be written as the *primal* LP\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "V(C) = & \\max \\hspace{2mm} v \\\\\n",
    "\\text{s.t. } \\hspace{2mm} v \\iota_n^T &\\leq x^T C  \\\\\n",
    "x &\\geq \\mathbf{0} \\\\\n",
    "\\iota_n^T x & = 1\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "while the *min-max problem* (the game from the minimizing player’s point\n",
    "of view) is the *dual* LP\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "V(C) = &\\min \\hspace{2mm} u \\\\\n",
    "\\text{s.t. } \\hspace{2mm}u \\iota_m &\\geq Cp  \\\\\n",
    "p &\\geq \\mathbf{0} \\\\\n",
    "\\iota_m^T p & = 1\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Hamburger, Thompson and Weil [[Hamburger *et al.*, 1967](https://python.quantecon.org/zreferences.html#id79)] view the input-output pair of the\n",
    "economy as payoff matrices of two-player zero-sum games.\n",
    "\n",
    "Using this interpretation, they restate Assumption I and II as follows\n",
    "\n",
    "$$\n",
    "V(-A) < 0\\quad\\quad \\text{and}\\quad\\quad V(B)>0\n",
    "$$\n",
    "\n",
    ">**Note**\n",
    ">\n",
    ">*Proof (Sketch)*:\n",
    "\n",
    "- $ \\Rightarrow $ $ V(B)>0 $ implies\n",
    "  $ x_0^T B \\gg \\mathbf{0} $, where $ x_0 $ is a maximizing\n",
    "  vector. Since $ B $ is non-negative, this requires that each\n",
    "  column of $ B $ has at least one positive entry, which is\n",
    "  Assumption I.  \n",
    "- $ \\Leftarrow $ From Assumption I and the fact\n",
    "  that $ p>\\mathbf{0} $, it follows that $ Bp > \\mathbf{0} $.\n",
    "  This implies that the maximizing player can always choose $ x $\n",
    "  so that $ x^TBp>0 $ so that it must be the case\n",
    "  that $ V(B)>0 $.  \n",
    "\n",
    "\n",
    "In order to (re)state Theorem I in terms of a particular two-player\n",
    "zero-sum game, we define a matrix for $ \\gamma\\in\\mathbb{R} $\n",
    "\n",
    "$$\n",
    "M(\\gamma) \\equiv B - \\gamma A\n",
    "$$\n",
    "\n",
    "For fixed $ \\gamma $, treating $ M(\\gamma) $ as a matrix game,\n",
    "calculating the solution of the game implies\n",
    "\n",
    "- If $ \\gamma > \\alpha_0 $, then for all $ x>0 $, there\n",
    "  $ \\exists j\\in\\{1, \\dots, n\\} $, s.t.\n",
    "  $ [x^T M(\\gamma)]_j < 0 $ implying\n",
    "  that $ V(M(\\gamma)) < 0 $.  \n",
    "- If $ \\gamma < \\beta_0 $, then for all $ p>0 $, there\n",
    "  $ \\exists i\\in\\{1, \\dots, m\\} $, s.t.\n",
    "  $ [M(\\gamma)p]_i > 0 $ implying that $ V(M(\\gamma)) > 0 $.  \n",
    "- If $ \\gamma \\in \\{\\beta_0, \\alpha_0\\} $, then (by Theorem I) the\n",
    "  optimal intensity and price vectors $ x_0 $ and $ p_0 $\n",
    "  satisfy  \n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "x_0^T M(\\gamma) \\geq \\mathbf{0}^T \\quad \\quad \\text{and}\\quad\\quad M(\\gamma) p_0 \\leq \\mathbf{0}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "That is, $ (x_0, p_0, 0) $ is a solution of the game\n",
    "$ M(\\gamma) $ so\n",
    "that $ V\\left(M(\\beta_0)\\right) = V\\left(M(\\alpha_0)\\right) = 0 $.\n",
    "\n",
    "- If $ \\beta_0 < \\alpha_0 $ and\n",
    "  $ \\gamma \\in (\\beta_0, \\alpha_0) $, then $ V(M(\\gamma)) = 0 $.  \n",
    "\n",
    "\n",
    "Moreover, if $ x' $ is optimal for the maximizing player in\n",
    "$ M(\\gamma') $ for $ \\gamma'\\in(\\beta_0, \\alpha_0) $ and\n",
    "$ p'' $ is optimal for the minimizing player in $ M(\\gamma'') $\n",
    "where $ \\gamma''\\in(\\beta_0, \\gamma') $, then $ (x', p'', 0) $\n",
    "is a solution for $ M(\\gamma) $ $ \\forall \\gamma\\in (\\gamma'', \\gamma') $.\n",
    "\n",
    "*Proof (Sketch):* If $ x' $ is optimal for a maximizing player in\n",
    "game $ M(\\gamma') $, then $ (x')^T M(\\gamma')\\geq \\mathbf{0}^T $ and so for all $ \\gamma<\\gamma' $.\n",
    "\n",
    "$$\n",
    "(x')^T M(\\gamma) = (x')^T M(\\gamma') + (x')^T(\\gamma' - \\gamma)A \\geq \\mathbf{0}^T\n",
    "$$\n",
    "\n",
    "hence $ V(M(\\gamma))\\geq 0 $. If $ p'' $ is optimal for a\n",
    "minimizing player in game $ M(\\gamma'') $, then $ M(\\gamma)p \\leq \\mathbf{0} $\n",
    "and so for all $ \\gamma''<\\gamma $\n",
    "\n",
    "$$\n",
    "M(\\gamma)p'' = M(\\gamma'') + (\\gamma'' - \\gamma)Ap'' \\leq \\mathbf{0}\n",
    "$$\n",
    "\n",
    "hence $ V(M(\\gamma))\\leq 0 $.\n",
    "\n",
    "It is clear from the above argument that $ \\beta_0 $, $ \\alpha_0 $ are the minimal and maximal $ \\gamma $ for which\n",
    "$ V(M(\\gamma))=0 $.\n",
    "\n",
    "Furthermore, Hamburger et al. [[Hamburger *et al.*, 1967](https://python.quantecon.org/zreferences.html#id79)] show that the\n",
    "function $ \\gamma \\mapsto V(M(\\gamma)) $ is continuous and\n",
    "nonincreasing in $ \\gamma $.\n",
    "\n",
    "This suggests an algorithm to compute\n",
    "$ (\\alpha_0, x_0) $ and $ (\\beta_0, p_0) $ for a given\n",
    "input-output pair $ (A, B) $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b6620ef",
   "metadata": {},
   "source": [
    "### Algorithm\n",
    "\n",
    "Hamburger, Thompson and Weil [[Hamburger *et al.*, 1967](https://python.quantecon.org/zreferences.html#id79)] propose a simple bisection algorithm\n",
    "to find the minimal and maximal roots (i.e. $ \\beta_0 $ and\n",
    "$ \\alpha_0 $) of the function $ \\gamma \\mapsto V(M(\\gamma)) $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dae4e2ff",
   "metadata": {},
   "source": [
    "#### Step 1\n",
    "\n",
    "First, notice that we can easily find trivial upper and lower bounds for\n",
    "$ \\alpha_0 $ and $ \\beta_0 $.\n",
    "\n",
    "- TEP requires that\n",
    "  $ x^T(B-\\alpha A)\\geq \\mathbf{0}^T $ and $ x > \\mathbf{0} $, so\n",
    "  if $ \\alpha $ is so large that\n",
    "  $ \\max_i\\{[(B-\\alpha A)\\iota_n]_i\\} < 0 $, then TEP ceases to have a\n",
    "  solution.  \n",
    "\n",
    "\n",
    "Accordingly, let **`UB`** be the $ \\alpha^{*} $ that\n",
    "solves $ \\max_i\\{[(B-\\alpha^{*} A)\\iota_n]_i\\} = 0 $.\n",
    "\n",
    "- Similar to\n",
    "  the upper bound, if $ \\beta $ is so low that\n",
    "  $ \\min_j\\{[\\iota^T_m(B-\\beta A)]_j\\}>0 $, then the EEP has no\n",
    "  solution and so we can define **`LB`** as the $ \\beta^{*} $ that\n",
    "  solves $ \\min_j\\{[\\iota^T_m(B-\\beta^{*} A)]_j\\}=0 $.  \n",
    "\n",
    "\n",
    "The *bounds* method calculates these trivial bounds for us"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca2cc46b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "n1.bounds()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "afb98ee2",
   "metadata": {},
   "source": [
    "#### Step 2\n",
    "\n",
    "Compute $ \\alpha_0 $ and $ \\beta_0 $\n",
    "\n",
    "- Finding $ \\alpha_0 $  \n",
    "  1. Fix $ \\gamma = \\frac{UB + LB}{2} $ and compute the solution\n",
    "    of the two-player zero-sum game associated\n",
    "    with $ M(\\gamma) $. We can use either the primal or the dual\n",
    "    LP problem.  \n",
    "  1. If $ V(M(\\gamma)) \\geq 0 $, then set $ LB = \\gamma $,\n",
    "    otherwise let $ UB = \\gamma $.  \n",
    "  1. Iterate on 1. and 2. until $ |UB - LB| < \\epsilon $.  \n",
    "- Finding $ \\beta_0 $  \n",
    "  1. Fix $ \\gamma = \\frac{UB + LB}{2} $ and compute the solution\n",
    "    of the two-player zero-sum game associated.\n",
    "    with $ M(\\gamma) $. We can use either the primal or the dual\n",
    "    LP problem.  \n",
    "  1. If $ V(M(\\gamma)) > 0 $, then set $ LB = \\gamma $,\n",
    "    otherwise let $ UB = \\gamma $.  \n",
    "  1. Iterate on 1. and 2. until $ |UB - LB| < \\epsilon $.  \n",
    "- *Existence*: Since $ V(M(LB))>0 $ and $ V(M(UB))<0 $ and\n",
    "  $ V(M(\\cdot)) $ is a continuous, nonincreasing function, there is\n",
    "  at least one $ \\gamma\\in[LB, UB] $, s.t. $ V(M(\\gamma))=0 $.  \n",
    "\n",
    "\n",
    "The *zerosum* method calculates the value and optimal strategies\n",
    "associated with a given $ \\gamma $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80d0a447",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "γ = 2\n",
    "\n",
    "print(f'Value of the game with γ = {γ}')\n",
    "print(n1.zerosum(γ=γ)[0])\n",
    "print('Intensity vector (from the primal)')\n",
    "print(n1.zerosum(γ=γ)[1])\n",
    "print('Price vector (from the dual)')\n",
    "print(n1.zerosum(γ=γ, dual=True)[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a500cf4",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "numb_grid = 100\n",
    "γ_grid = np.linspace(0.4, 2.1, numb_grid)\n",
    "\n",
    "value_ex1_grid = np.asarray([n1.zerosum(γ=γ_grid[i])[0]\n",
    "                            for i in range(numb_grid)])\n",
    "value_ex2_grid = np.asarray([n2.zerosum(γ=γ_grid[i])[0]\n",
    "                            for i in range(numb_grid)])\n",
    "\n",
    "fig, axes = plt.subplots(1, 2, figsize=(14, 5), sharey=True)\n",
    "fig.suptitle(r'The function $V(M(\\gamma))$', fontsize=16)\n",
    "\n",
    "for ax, grid, N, i in zip(axes, (value_ex1_grid, value_ex2_grid),\n",
    "                          (n1, n2), (1, 2)):\n",
    "    ax.plot(γ_grid, grid)\n",
    "    ax.set(title=f'Example {i}', xlabel='$\\gamma$')\n",
    "    ax.axhline(0, c='k', lw=1)\n",
    "    ax.axvline(N.bounds()[0], c='r', ls='--', label='lower bound')\n",
    "    ax.axvline(N.bounds()[1], c='g', ls='--', label='upper bound')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9612417",
   "metadata": {},
   "source": [
    "The *expansion* method implements the bisection algorithm for\n",
    "$ \\alpha_0 $ (and uses the primal LP problem for $ x_0 $)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c85dd21a",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "α_0, x, p = n1.expansion()\n",
    "print(f'α_0 = {α_0}')\n",
    "print(f'x_0 = {x}')\n",
    "print(f'The corresponding p from the dual = {p}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08d1d955",
   "metadata": {},
   "source": [
    "The *interest* method implements the bisection algorithm for\n",
    "$ \\beta_0 $ (and uses the dual LP problem for $ p_0 $)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d77ef1bb",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "β_0, x, p = n1.interest()\n",
    "print(f'β_0 = {β_0}')\n",
    "print(f'p_0 = {p}')\n",
    "print(f'The corresponding x from the primal = {x}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bc8174f",
   "metadata": {},
   "source": [
    "Of course, when $ \\gamma^* $ is unique, it is irrelevant which one of the two methods we use  – both work.\n",
    "\n",
    "In particular, as will be shown below, in\n",
    "case of an irreducible $ (A,B) $ (like in Example 1), the maximal\n",
    "and minimal roots of $ V(M(\\gamma)) $ necessarily coincide implying\n",
    "a ‘‘full duality’’ result, i.e. $ \\alpha_0 = \\beta_0 = \\gamma^* $ so that the expansion (and interest) rate $ \\gamma^* $ is unique."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8035f07b",
   "metadata": {},
   "source": [
    "### Uniqueness and Irreducibility\n",
    "\n",
    "As an illustration, compute first the maximal and minimal roots of\n",
    "$ V(M(\\cdot)) $ for our Example 2 that has a reducible\n",
    "input-output pair $ (A, B) $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de43b849",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "α_0, x, p = n2.expansion()\n",
    "print(f'α_0 = {α_0}')\n",
    "print(f'x_0 = {x}')\n",
    "print(f'The corresponding p from the dual = {p}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac4be00d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "β_0, x, p = n2.interest()\n",
    "print(f'β_0 = {β_0}')\n",
    "print(f'p_0 = {p}')\n",
    "print(f'The corresponding x from the primal = {x}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bddcfc8",
   "metadata": {},
   "source": [
    "As we can see, with a reducible $ (A,B) $, the roots found by the\n",
    "bisection algorithms might differ, so there might be multiple\n",
    "$ \\gamma^* $ that make the value of the game\n",
    "with $ M(\\gamma^*) $ zero. (see the figure above).\n",
    "\n",
    "Indeed, although the von Neumann theorem assures existence of the\n",
    "equilibrium, Assumptions I and II are not sufficient for uniqueness.\n",
    "Nonetheless, Kemeny et al. (1967) show that there are at most finitely\n",
    "many economic solutions, meaning that there are only finitely many\n",
    "$ \\gamma^* $ that satisfy $ V(M(\\gamma^*)) = 0 $ and\n",
    "$ x_0^TBp_0 > 0 $ and that for each such $ \\gamma^*_i $, there\n",
    "is a self-contained part of the economy (a sub-economy) that in\n",
    "equilibrium can expand independently with the expansion\n",
    "coefficient $ \\gamma^*_i $.\n",
    "\n",
    "The following theorem (see Theorem 9.10. in Gale [[Gale, 1989](https://python.quantecon.org/zreferences.html#id80)]) asserts that\n",
    "imposing irreducibility is sufficient for uniqueness of\n",
    "$ (\\gamma^*, x_0, p_0) $.\n",
    "\n",
    "**Theorem II:** Adopt the conditions of Theorem 1. If the economy\n",
    "$ (A,B) $ is irreducible, then $ \\gamma^*=\\alpha_0=\\beta_0 $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eaa49079",
   "metadata": {},
   "source": [
    "### A Special Case\n",
    "\n",
    "There is a special $ (A,B) $ that allows us to simplify the solution\n",
    "method significantly by invoking the powerful Perron-Frobenius theorem\n",
    "for non-negative matrices.\n",
    "\n",
    "**Definition:** We call an economy *simple* if it satisfies\n",
    "\n",
    "- $ n=m $  \n",
    "- Each activity produces exactly one good  \n",
    "- Each good is produced by one and only one activity.  \n",
    "\n",
    "\n",
    "These assumptions imply that $ B=I_n $, i.e., that $ B $ can be\n",
    "written as an identity matrix (possibly after reshuffling its rows and\n",
    "columns).\n",
    "\n",
    "The simple model has the following special property (Theorem 9.11. in Gale [[Gale, 1989](https://python.quantecon.org/zreferences.html#id80)]): if $ x_0 $ and $ \\alpha_0>0 $ solve the TEP\n",
    "with $ (A,I_n) $, then\n",
    "\n",
    "$$\n",
    "x_0^T = \\alpha_0 x_0^T A\\hspace{1cm}\\Leftrightarrow\\hspace{1cm}x_0^T\n",
    "A=\\left(\\frac{1}{\\alpha_0}\\right)x_0^T\n",
    "$$\n",
    "\n",
    "The latter shows that $ 1/\\alpha_0 $ is a positive eigenvalue of\n",
    "$ A $ and $ x_0 $ is the corresponding non-negative left\n",
    "eigenvector.\n",
    "\n",
    "The classic result of **Perron and Frobenius** implies\n",
    "that a non-negative matrix  has a non-negative\n",
    "eigenvalue-eigenvector pair.\n",
    "\n",
    "Moreover, if $ A $ is irreducible, then\n",
    "the optimal intensity vector $ x_0 $ is positive and *unique* up to\n",
    "multiplication by a positive scalar.\n",
    "\n",
    "Suppose that $ A $ is reducible with $ k $ irreducible subsets\n",
    "$ S_1,\\dots,S_k $. Let $ A_i $ be the submatrix corresponding to\n",
    "$ S_i $ and let $ \\alpha_i $ and $ \\beta_i $ be the\n",
    "associated expansion and interest factors, respectively. Then we have\n",
    "\n",
    "$$\n",
    "\\alpha_0 = \\max_i \\{\\alpha_i\\}\\hspace{1cm}\\text{and}\\hspace{1cm}\\beta_0 = \\min_i \\{\\beta_i\\}\n",
    "$$"
   ]
  }
 ],
 "metadata": {
  "date": 1714442510.9304576,
  "filename": "von_neumann_model.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Von Neumann Growth Model (and a Generalization)"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}