BCG_incomplete_mkts.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a8ede616",
   "metadata": {},
   "source": [
    "\n",
    "<a id='bcg-incomplete-final'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81fd163f",
   "metadata": {},
   "source": [
    "# Equilibrium Capital Structures with Incomplete Markets\n",
    "\n",
    "In addition to what’s in Anaconda, this lecture will need the following libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57879a5f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "!pip install --upgrade quantecon\n",
    "!conda install -y -c plotly plotly plotly-orca"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9713d1da",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "This is an extension of an earlier lecture [Irrelevance of Capital Structure with Complete Markets](https://python-advanced.quantecon.org/BCG_complete_mkts.html) about a **complete markets**\n",
    "model.\n",
    "\n",
    "In contrast to that lecture, this one  describes an instance of a model authored by Bisin, Clementi, and Gottardi [[Bisin *et al.*, 2018](https://python-advanced.quantecon.org/zreferences.html#id32)]\n",
    "in which financial markets are **incomplete**.\n",
    "\n",
    "Instead of being able to trade equities and a full set of one-period\n",
    "Arrow securities as they can in [Irrelevance of Capital Structure with Complete Markets](https://python-advanced.quantecon.org/BCG_complete_mkts.html), here consumers and firms trade only equity and a bond.\n",
    "\n",
    "It is useful to watch how outcomes differ  in the two settings.\n",
    "\n",
    "In the complete markets economy in [Irrelevance of Capital Structure with Complete Markets](https://python-advanced.quantecon.org/BCG_complete_mkts.html)\n",
    "\n",
    "- there is a unique  stochastic discount factor that prices all assets  \n",
    "- consumers’ portfolio choices are indeterminate  \n",
    "- firms’ financial structures are indeterminate, so the model embodies an instance of a Modigliani-Miller irrelevance theorem [[Modigliani and Miller, 1958](https://python-advanced.quantecon.org/zreferences.html#id33)]  \n",
    "- the aggregate of all firms’ financial structures are indeterminate, a consequence of there being redundant assets  \n",
    "\n",
    "\n",
    "In the incomplete markets economy studied here\n",
    "\n",
    "- there is a not a unique equilibrium stochastic discount factor  \n",
    "- different stochastic discount factors price different assets  \n",
    "- consumers’ portfolio choices are determinate  \n",
    "- while **individual** firms’ financial structures are indeterminate, thus conforming to part of a Modigliani-Miller theorem,\n",
    "  [[Modigliani and Miller, 1958](https://python-advanced.quantecon.org/zreferences.html#id33)], the **aggregate**  of all firms’ financial structures **is** determinate.  \n",
    "\n",
    "\n",
    "A  `Big K, little k` analysis played an important role in the  previous lecture [Irrelevance of Capital Structure with Complete Markets](https://python-advanced.quantecon.org/BCG_complete_mkts.html).\n",
    "\n",
    "A more subtle version of a `Big K, little k`   features in  the   BCG incomplete markets environment here.\n",
    "\n",
    "We use it to convey the heart of what BCG call a **rational conjectures** equilibrium in which conjectures are about\n",
    "equilibrium pricing functions in regions of the state space that  an average consumer or firm does not visit in equilibrium.\n",
    "\n",
    "Note that the absence of complete markets means that now we cannot compute competitive equilibrium prices and allocations by first solving the simple planning problem that we did in [Irrelevance of Capital Structure with Complete Markets](https://python-advanced.quantecon.org/BCG_complete_mkts.html).\n",
    "\n",
    "Instead, we compute an equilibrium by solving a system of simultaneous inequalities.\n",
    "\n",
    "(Here we do not address the interesting question of whether there is a *different* planning problem that we could use to compute a\n",
    "competitive equlibrium allocation.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c5905d6",
   "metadata": {},
   "source": [
    "### Setup\n",
    "\n",
    "We adopt specifications of preferences and technologies used by Bisin,\n",
    "Clemente, and Gottardi (2018) [[Bisin *et al.*, 2018](https://python-advanced.quantecon.org/zreferences.html#id32)] and in our earlier lecture on a complete markets\n",
    "version of their model.\n",
    "\n",
    "The economy lasts for two periods, $ t=0, 1 $.\n",
    "\n",
    "There are two types of consumers named $ i=1,2 $.\n",
    "\n",
    "A scalar random variable $ \\epsilon $ affects both\n",
    "\n",
    "- a representative firm’s physical return $ f(k)e^\\epsilon $ in\n",
    "  period $ 1 $ from investing $ k \\geq 0 $ in capital in period\n",
    "  $ 0 $.  \n",
    "- period $ 1 $ endowments $ w_1^i(\\epsilon) $ of the\n",
    "  consumption good for agents $ i =1 $ and $ i=2 $.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "669b916a",
   "metadata": {},
   "source": [
    "### Ownership\n",
    "\n",
    "A consumer of type $ i $ is endowed with $ w_0^i $ units of the\n",
    "time $ 0 $ good and $ w_1^i(\\epsilon) $ of the time $ 1 $\n",
    "good when the random variable takes value $ \\epsilon $.\n",
    "\n",
    "At the start of period $ 0 $, a consumer of type $ i $ also owns\n",
    "$ \\theta^i_0 $ shares of a representative firm."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "362957da",
   "metadata": {},
   "source": [
    "### Measures of agents and firms\n",
    "\n",
    "As in the companion lecture [Irrelevance of Capital Structure with Complete Markets](https://python-advanced.quantecon.org/BCG_complete_mkts.html) that studies a complete markets version of\n",
    "the model, we follow BCG in assuming that there are unit measures of\n",
    "\n",
    "- consumers of type $ i=1 $  \n",
    "- consumers of type $ i=2 $  \n",
    "- firms with access to a production technology that converts\n",
    "  $ k $ units of time $ 0 $ good into\n",
    "  $ A k^\\alpha e^\\epsilon $ units of the time $ 1 $ good in\n",
    "  random state $ \\epsilon $  \n",
    "\n",
    "\n",
    "Thus, let $ \\omega \\in [0,1] $ index a particular consumer of type\n",
    "$ i $.\n",
    "\n",
    "Then define Big $ C^i $ as\n",
    "\n",
    "$$\n",
    "C^i = \\int_0^1 c^i(\\omega) d \\, \\omega\n",
    "$$\n",
    "\n",
    "with components\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "C^i_0  & =  \\int_0^1 c^i_0(\\omega) d \\, \\omega \\cr\n",
    "C^i_1(\\epsilon) &  = \\int_0^1 c^i_1(\\epsilon;\\omega) d \\, \\omega\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "In the same spirit, let $ \\zeta \\in [0,1] $ index a particular firm\n",
    "and let firm $ \\zeta $ purchase $ k(\\zeta) $ units of capital\n",
    "and issue $ b(\\zeta) $ bonds.\n",
    "\n",
    "Then define Big $ K $ and Big $ B $ as\n",
    "\n",
    "$$\n",
    "K = \\int_0^1 k(\\zeta) d \\, \\zeta, \\quad B = \\int_0^1 b(\\zeta) d \\, \\zeta\n",
    "$$\n",
    "\n",
    "The assumption that there are equal measures of our three types of\n",
    "agents justifies our assumption that  each individual agent is a\n",
    "powerless **price taker**:\n",
    "\n",
    "- an individual consumer chooses its own (infinitesimal) part\n",
    "  $ c^i(\\omega) $ of $ C^i $ taking prices as given  \n",
    "- an individual firm chooses its own (infinitesmimal) part\n",
    "  $ k(\\zeta) $ of $ K $ and $ b(\\zeta) $ of $ B $\n",
    "  taking pricing functions as given  \n",
    "- However, equilibrium prices depend on the `Big K, Big B, Big C`\n",
    "  objects $ K $, $ B $, and $ C $  \n",
    "\n",
    "\n",
    "The assumption about measures of agents is a powerful device for making\n",
    "a host of competitive agents take as given the equilibrium prices that\n",
    "turn out to be determined by the decisions of hosts of agents who are just like\n",
    "them.\n",
    "\n",
    "We call an equilibrium **symmetric** if\n",
    "\n",
    "- all type $ i $ consumers choose the same consumption profiles so\n",
    "  that $ c^i(\\omega) = C^i $ for all $ \\omega \\in [0,1] $  \n",
    "- all firms choose the same levels of $ k $ and $ b $ so that\n",
    "  $ k(\\zeta) = K $, $ b(\\zeta) = B $ for all\n",
    "  $ \\zeta \\in [0,1] $  \n",
    "\n",
    "\n",
    "In this lecture, we restrict ourselves to describing symmetric\n",
    "equilibria."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "102bfe36",
   "metadata": {},
   "source": [
    "### Endowments\n",
    "\n",
    "Per capital economy-wide endowments in periods $ 0 $ and $ 1 $ are\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "w_0 & = w_0^1 + w_0^2 \\cr\n",
    "w_1(\\epsilon) & =  w_1^1(\\epsilon) + w_1^2(\\epsilon) \\textrm{ in state }\\epsilon\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0d5dc7f",
   "metadata": {},
   "source": [
    "### Feasibility:\n",
    "\n",
    "Where $ \\alpha \\in (0,1) $ and $ A >0 $\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    " C_0^1 + C_0^2 & = w_0^1 + w_0^2 - K \\cr\n",
    " C_1^1(\\epsilon) + C_1^2(\\epsilon) & =  w_1^1(\\epsilon) + w_1^2(\\epsilon) + e^\\epsilon \\int_0^1 f(k(\\zeta)) d \\zeta, \\quad k \\geq 0\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "where $ f(k) = A k^\\alpha, A >0, \\alpha \\in (0,1) $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d656ffee",
   "metadata": {},
   "source": [
    "### Parameterizations\n",
    "\n",
    "Following BCG, we shall employ the following parameterizations:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\epsilon & \\sim {\\mathcal N}(\\mu, \\sigma^2) \\cr\n",
    "u(c) & = \\frac{c^{1-\\gamma}}{1 - \\gamma} \\cr\n",
    "w_1^i(\\epsilon) & = e^{- \\chi_i \\mu - .5 \\chi_i^2 \\sigma^2 + \\chi_i \\epsilon} , \\quad \\chi_i \\in [0,1]\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Sometimes instead of asuming $ \\epsilon \\sim g(\\epsilon) = {\\mathcal N}(0,\\sigma^2) $,\n",
    "we’ll assume that $ g(\\cdot) $ is a probability\n",
    "mass function that serves as a discrete approximation to a standardized\n",
    "normal density."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b34c8f3",
   "metadata": {},
   "source": [
    "### Preferences:\n",
    "\n",
    "A consumer of type $ i $ orders period $ 0 $ consumption\n",
    "$ c_0^i $ and state $ \\epsilon $-period $ 1 $ consumption\n",
    "$ c^i(\\epsilon) $ by\n",
    "\n",
    "$$\n",
    "u^i = u(c_0^i) + \\beta \\int u(c_1^i(\\epsilon)) g (\\epsilon) d \\epsilon, \\quad i = 1,2\n",
    "$$\n",
    "\n",
    "$ \\beta \\in (0,1) $ and the one-period utility function is\n",
    "\n",
    "$$\n",
    "u(c) = \\begin{cases}\n",
    "\\frac{c^{1 -\\gamma}} { 1 - \\gamma} & \\textrm{if  } \\gamma \\neq 1 \\\\\n",
    "\\log c & \\textrm{if  } \\gamma = 1\n",
    "\\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b32affe",
   "metadata": {},
   "source": [
    "### Risk-sharing motives\n",
    "\n",
    "The two types of agents’ period $ 1 $ endowments have different correlations with\n",
    "the physical return on capital.\n",
    "\n",
    "Endowment differences give agents incentives to trade risks that in the\n",
    "complete market version of the model showed up in their demands for\n",
    "equity and in their demands and supplies of one-period Arrow securities.\n",
    "\n",
    "In the incomplete-markets setting under study here, these differences\n",
    "show up in differences in the two types of consumers’ demands for a\n",
    "typical firm’s bonds and equity, the only two assets that agents can now\n",
    "trade."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b684f66",
   "metadata": {},
   "source": [
    "## Asset Markets\n",
    "\n",
    "Markets are incomplete: *ex cathedra* we the model builders declare that only equities and bonds issued by representative\n",
    "firms can be traded.\n",
    "\n",
    "Let $ \\theta^i $ and $ \\xi^i $ be a consumer of type\n",
    "$ i $’s post-trade holdings of equity and bonds, respectively.\n",
    "\n",
    "A firm issues bonds promising to pay $ b $ units of consumption at\n",
    "time $ t=1 $ and purchases $ k $ units of physical capital at\n",
    "time $ t=0 $.\n",
    "\n",
    "When $ e^\\epsilon A k^\\alpha < b $ at time $ 1 $, the firm defaults and its output is\n",
    "divided equally among bondholders.\n",
    "\n",
    "Evidently, when the productivity shock\n",
    "$ \\epsilon < \\epsilon^* = \\log \\left(\\frac{b}{ Ak^\\alpha}\\right) $,\n",
    "the firm defaults on its debt\n",
    "\n",
    "Payoffs to equity and debt at date 1 as functions of the productivity\n",
    "shock $ \\epsilon $ are thus\n",
    "\n",
    "\n",
    "<a id='equation-payofffns'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "d^e(k,b;\\epsilon) &= \\max \\left\\{ e^\\epsilon  A k^\\alpha - b, 0 \\right\\} \\\\\n",
    "d^b(k,b;\\epsilon) &= \\min \\left\\{ \\frac{e^\\epsilon  A k^\\alpha}{b}, 1 \\right\\}\n",
    "\\end{aligned} \\tag{39.1}\n",
    "$$\n",
    "\n",
    "A firm faces a bond price function $ p(k,b) $ when it issues\n",
    "$ b $ bonds and purchases $ k $ units of physical capital.\n",
    "\n",
    "A firm’s equity is worth $ q(k,b) $ when it issues $ b $ bonds\n",
    "and purchases $ k $ units of physical capital.\n",
    "\n",
    "A firm regards an equity-pricing function $ q(k,b) $ and a bond\n",
    "pricing function $ p(k,b) $ as exogenous in the sense that they are\n",
    "not affected by its choices of $ k $ and $ b $.\n",
    "\n",
    "Consumers face equilibrium prices $ \\check q $ and $ \\check p $\n",
    "for bonds and equities, where $ \\check q $ and $ \\check p $ are\n",
    "both scalars.\n",
    "\n",
    "Consumers are price takers and only need to know the scalars $ \\check q, \\check p $.\n",
    "\n",
    "Firms are *price function* takers and  must know the functions $ q(k,b), p(k,b) $ in order\n",
    "completely to pose their optimum problems."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0a63ef6",
   "metadata": {},
   "source": [
    "### Consumers\n",
    "\n",
    "Each consumer of type $ i $ is endowed with $ w_0^i $ of the\n",
    "time $ 0 $ consumption good, $ w_1^i(\\epsilon) $ of the time\n",
    "$ 1 $, state $ \\epsilon $ consumption good  and also owns a fraction\n",
    "$ \\theta^i_0 \\in (0,1) $ of the initial value of a representative\n",
    "firm, where $ \\theta^1_0 + \\theta^2_0 = 1 $.\n",
    "\n",
    "The initial value of a representative firm is $ V $ (an object to be\n",
    "determined in a rational expectations equilibrium).\n",
    "\n",
    "Consumer $ i $ buys $ \\theta^i $ shares of equity and buys bonds\n",
    "worth $ \\check p \\xi^i $ where $ \\check p $ is the bond price.\n",
    "\n",
    "Being a price-taker, a consumer takes $ V $, $ \\check q $, $ \\check p $, and $ K, B $\n",
    "as given.\n",
    "\n",
    "Consumers know that equilibrium payoff functions for bonds and equities take the form\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "d^e(K,B;\\epsilon) &= \\max \\left\\{ e^\\epsilon  A K^\\alpha - B, 0 \\right\\} \\\\\n",
    "d^b(K,B;\\epsilon) &= \\min \\left\\{ \\frac{e^\\epsilon  A K^\\alpha}{B}, 1 \\right\\}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Consumer $ i $’s optimization problem is\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\max_{c^i_0,\\theta^i,\\xi^i,c^i_1(\\epsilon)} & u(c^i_0) + \\beta \\int u(c^i(\\epsilon)) g(\\epsilon) \\ d\\epsilon \\\\\n",
    "\\text{subject to } \\quad\n",
    "& c^i_0 = w^i_0 + \\theta^i_0V - \\check q\\theta^i - \\check p \\xi^i, \\\\\n",
    "& c^i_1(\\epsilon) = w^i_1(\\epsilon) + \\theta^i d^e(K,B;\\epsilon) + \\xi^i d^b(K,B;\\epsilon) \\ \\forall \\ \\epsilon, \\\\\n",
    "& \\theta^i \\geq 0, \\xi^i \\geq 0.\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The last two inequalities impose that the consumer cannot short sell either\n",
    "equity or bonds.\n",
    "\n",
    "In a rational expectations equilibrium, $ \\check q = q(K,B) $ and $ \\check p = p(K,B) $\n",
    "\n",
    "We form consumer $ i $’s Lagrangian:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "L^i := &  u(c^i_0) + \\beta \\int u(c^i(\\epsilon)) g(\\epsilon) \\ d\\epsilon  \\\\\n",
    "     & +\\lambda^i_0 [w^i_0 + \\theta_0V - \\check q\\theta^i - \\check p \\xi^i - c^i_0]  \\\\\n",
    "     & +  \\beta \\int \\lambda^i_1(\\epsilon) \\left[ w^i_1(\\epsilon) + \\theta^i d^e(K,B;\\epsilon) + \\xi^i d^b(K,B;\\epsilon) - c^i_1(\\epsilon) \\right] g(\\epsilon) \\ d\\epsilon\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Consumer $ i $’s first-order necessary conditions for an optimum\n",
    "include:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "c^i_0:& \\quad u^\\prime(c^i_0) = \\lambda^i_0 \\\\\n",
    "c^i_1(\\epsilon):& \\quad u^\\prime(c^i_1(\\epsilon)) = \\lambda^i_1(\\epsilon) \\\\\n",
    "\\theta^i:& \\quad \\beta \\int \\lambda^i_1(\\epsilon) d^e(K,B;\\epsilon) g(\\epsilon) \\ d\\epsilon \\leq \\lambda^i_0 \\check q \\quad (= \\ \\ \\text{if} \\ \\ \\theta^i>0) \\\\\n",
    "\\xi^i:& \\quad \\beta \\int \\lambda^i_1(\\epsilon) d^b(K,B;\\epsilon) g(\\epsilon) \\ d\\epsilon \\leq \\lambda^i_0 \\check p \\quad (= \\ \\ \\text{if} \\ \\ b^i>0) \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "We can combine and rearrange consumer $ i $’s first-order\n",
    "conditions to become:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\check q \\geq \\beta \\int \\frac{u^\\prime(c^i_1(\\epsilon))}{u^\\prime(c^i_0)} d^e(K,B;\\epsilon) g(\\epsilon) \\ d\\epsilon \\quad (= \\ \\ \\text{if} \\ \\ \\theta^i>0) \\\\\n",
    "\\check p \\geq \\beta \\int \\frac{u^\\prime(c^i_1(\\epsilon))}{u^\\prime(c^i_0)} d^b(K,B;\\epsilon) g(\\epsilon) \\ d\\epsilon \\quad (= \\ \\ \\text{if} \\ \\ b^i>0)\\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "These inequalities imply that in a symmetric rational expectations equilibrium  consumption allocations and\n",
    "prices  satisfy\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\check q = \\max_i \\beta \\int \\frac{u^\\prime(c^i_1(\\epsilon))}{u^\\prime(c^i_0)} d^e(K,B;\\epsilon) g(\\epsilon) \\ d\\epsilon \\\\\n",
    "\\check p = \\max_i \\beta \\int \\frac{u^\\prime(c^i_1(\\epsilon))}{u^\\prime(c^i_0)} d^b(K,B;\\epsilon) g(\\epsilon) \\ d\\epsilon \\\\\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2008e2ce",
   "metadata": {},
   "source": [
    "### Pricing functions\n",
    "\n",
    "When individual firms solve their optimization problems, they take big\n",
    "$ C^i $’s as fixed objects that they don’t influence.\n",
    "\n",
    "A representative firm faces a price function $ q(k,b) $ for its\n",
    "equity and a price function $ p(k, b) $ per unit of bonds that\n",
    "satisfy\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "q(k,b) = \\max_i \\beta \\int \\frac{u^\\prime(C^i_1(\\epsilon))}{u^\\prime(C^i_0)} d^e(k,b;\\epsilon) g(\\epsilon) \\ d\\epsilon \\\\\n",
    "p(k,b) = \\max_i \\beta \\int \\frac{u^\\prime(C^i_1(\\epsilon))}{u^\\prime(C^i_0)} d^b(k,b;\\epsilon) g(\\epsilon) \\ d\\epsilon \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "where the payoff functions are described by equations [(39.1)](#equation-payofffns).\n",
    "\n",
    "Notice the appearance of big $ C^i $’s on the right sides of these\n",
    "two equations that define equilibrium pricing functions.\n",
    "\n",
    "The two price functions describe outcomes not only for equilibrium choices\n",
    "$ K, B $ of capital $ k $ and debt $ b $, but also for any\n",
    "**out-of-equilibrium** pairs $ (k, b) \\neq (K, B) $.\n",
    "\n",
    "The firm is assumed to know both price functions.\n",
    "\n",
    "This means that the firm understands that its choice of $ k,b $ influences how markets price its equity and debt.\n",
    "\n",
    "This package of assumptions is sometimes called  **rational conjectures** (about price functions).\n",
    "\n",
    "BCG give credit to Makowski for emphasizing and clarifying how rational conjectures are components of  rational expectations equilibria."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71df309e",
   "metadata": {},
   "source": [
    "### Firms\n",
    "\n",
    "The firm chooses capital $ k $ and debt $ b $ to maximize its\n",
    "market value:\n",
    "\n",
    "$$\n",
    "V \\equiv \\max_{k,b} -k + q(k,b) + p(k,b) b\n",
    "$$\n",
    "\n",
    "Attributing value maximization to the firm is a good idea because in equilibrium consumers of both types\n",
    "*want* a firm to maximize its value.\n",
    "\n",
    "In the special quantitative examples studied here\n",
    "\n",
    "- consumers of  types $ i=1,2 $ both  hold equity  \n",
    "- only consumers of type $ i=2 $ hold debt; consumers of type\n",
    "  $ i=1 $ hold none.  \n",
    "\n",
    "\n",
    "These outcomes occur because we follow BCG and set parameters so that a\n",
    "type 2 consumer’s stochastic endowment of the consumption good in period\n",
    "$ 1 $ is more correlated with the firm’s output than is a type 1\n",
    "consumer’s.\n",
    "\n",
    "This gives consumers of type $ 2 $ a motive to hedge their second period\n",
    "endowment risk by holding bonds (they also choose to\n",
    "hold some equity).\n",
    "\n",
    "These outcomes mean that the pricing functions end up\n",
    "satisfying\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "q(k,b) &= \\beta \\int \\frac{u^\\prime(C^1_1(\\epsilon))}{u^\\prime(C^1_0)} d^e(k,b;\\epsilon) g(\\epsilon) \\ d\\epsilon = \\beta \\int \\frac{u^\\prime(C^2_1(\\epsilon))}{u^\\prime(C^2_0)} d^e(k,b;\\epsilon) g(\\epsilon) \\ d\\epsilon \\\\\n",
    "p(k,b) &= \\beta \\int \\frac{u^\\prime(C^2_1(\\epsilon))}{u^\\prime(C^2_0)} d^b(k,b;\\epsilon) g(\\epsilon) \\ d\\epsilon \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Recall that\n",
    "$ \\epsilon^*(k,b) \\equiv \\log\\left(\\frac{b}{Ak^\\alpha}\\right) $ is a\n",
    "firm’s default threshold.\n",
    "\n",
    "We can rewrite the pricing functions as:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "q(k,b) &= \\beta \\int_{\\epsilon^*}^\\infty \\frac{u^\\prime(C^i_1(\\epsilon))}{u^\\prime(C^i_0)} \\left( e^\\epsilon Ak^\\alpha - b \\right) g(\\epsilon) \\ d\\epsilon, \\quad i=1,2\\\\\n",
    "p(k,b) &= \\beta \\int^{\\epsilon^*}_{-\\infty} \\frac{u^\\prime(C^2_1(\\epsilon))}{u^\\prime(C^2_0)} \\left( \\frac{e^\\epsilon Ak^\\alpha}{b} \\right) g(\\epsilon) \\ d\\epsilon + \\beta \\int_{\\epsilon^*}^{\\infty} \\frac{u^\\prime(C^2_1(\\epsilon))}{u^\\prime(C^2_0)} g(\\epsilon) \\ d\\epsilon \\\\\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b50774ad",
   "metadata": {},
   "source": [
    "#### Firm’s optimization problem\n",
    "\n",
    "The firm’s optimization problem is\n",
    "\n",
    "$$\n",
    "V \\equiv \\max_{k,b} \\left\\{ -k + q(k,b) + p(k, b) b \\right\\}\n",
    "$$\n",
    "\n",
    "The firm’s first-order necessary conditions with respect to $ k $\n",
    "and $ b $, respectively, are\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "k: \\quad &   -1 + \\frac{\\partial q(k,b)}{\\partial k} + b \\frac{\\partial p(q,b)}{\\partial k}  = 0 \\cr\n",
    "    b: \\quad &  \\frac{\\partial q(k,b)}{\\partial b} + p(k,b) + b \\frac{\\partial p(k,b)}{\\partial b} = 0\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "We use the Leibniz integral rule several times to arrive at\n",
    "the following derivatives:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial q(k,b)}{\\partial k}  = \\beta \\alpha A k^{\\alpha-1} \\int_{\\epsilon^*}^\\infty \\frac{u'(C_1^i(\\epsilon))}{u'(C_0^i)}\n",
    "     e^\\epsilon g(\\epsilon) d \\epsilon, \\quad i=1,2\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\frac{\\partial q(k,b)}{\\partial b}  = -\\beta \\int_{\\epsilon^*}^\\infty \\frac{u'(C_1^i(\\epsilon))}{u'(C_0^i)} g(\\epsilon) d \\epsilon, \\quad i=1,2\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\frac{\\partial p(k,b)}{\\partial k} =  \\beta \\alpha \\frac{A k^{\\alpha -1}}{b} \\int_{-\\infty}^{\\epsilon^*} \\frac{u'(C_1^2(\\epsilon))}{u'(C_0^2)}   g(\\epsilon) d \\epsilon\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\frac{\\partial p(k,b)}{\\partial b} = - \\beta \\frac{A k^\\alpha}{b^2} \\int_{-\\infty}^{\\epsilon^*}  \\frac{u'(C_1^2(\\epsilon))}{u'(C_0^2)} e^\\epsilon  g(\\epsilon) d \\epsilon\n",
    "$$\n",
    "\n",
    "**Special case:** We confine ourselves to a special case in which both types of\n",
    "consumer hold positive equities so that\n",
    "$ \\frac{\\partial q(k,b)}{\\partial k} $ and\n",
    "$ \\frac{\\partial q(k,b)}{\\partial b} $ are related to rates of\n",
    "intertemporal substitution for both agents.\n",
    "\n",
    "Substituting these partial derivatives into the above first-order\n",
    "conditions for $ k $ and $ b $, respectively, we obtain the\n",
    "following versions of those first order conditions:\n",
    "\n",
    "\n",
    "<a id='equation-eqn1'></a>\n",
    "$$\n",
    "k: \\quad -1 + \\beta \\alpha A k^{\\alpha -1} \\int_{-\\infty}^\\infty \\frac{u'(C_1^2(\\epsilon))}{u'(C_0^2)} e^\\epsilon  g(\\epsilon) d \\epsilon = 0 \\tag{39.2}\n",
    "$$\n",
    "\n",
    "\n",
    "<a id='equation-eqn2'></a>\n",
    "$$\n",
    "b:  \\quad\n",
    "\\int_{\\epsilon^*}^\\infty \\left( \\frac{u^\\prime(C^1_1(\\epsilon))}{u^\\prime(C^1_0)} \\right) g(\\epsilon) \\ d\\epsilon = \\int_{\\epsilon^*}^\\infty \\left( \\frac{u^\\prime(C^2_1(\\epsilon))}{u^\\prime(C^2_0)} \\right)  g(\\epsilon) \\ d\\epsilon \\tag{39.3}\n",
    "$$\n",
    "\n",
    "where again recall that\n",
    "$ \\epsilon^*(k,b) \\equiv \\log\\left(\\frac{b}{Ak^\\alpha}\\right) $.\n",
    "\n",
    "Taking $ C_0^i, C_1^i(\\epsilon) $ as given, these are two equations\n",
    "that we want to solve for the firm’s optimal decisions $ k, b $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df250c19",
   "metadata": {},
   "source": [
    "## Equilibrium verification\n",
    "\n",
    "On page 5 of BCG (2018), the authors say\n",
    "\n",
    "*If the price conjectures corresponding to the plan chosen by firms in\n",
    "equilibrium are correct, that is equal to the market prices* $ \\check q $ *and* $ \\check p $, *it is immediate to verify that\n",
    "the rationality of the conjecture coincides with the agents’ Euler\n",
    "equations.*\n",
    "\n",
    "Here BCG are describing how they go about verifying that when they set\n",
    "little $ k $, little $ b $ from the firm’s first-order\n",
    "conditions equal to the big $ K $, big $ B $ at the big\n",
    "$ C $’s that appear in the pricing functions, then\n",
    "\n",
    "- consumers’ Euler equations are satisfied if little $ c $’s are\n",
    "  equated to  Big $ C $’s  \n",
    "- firms’ first-order necessary conditions for $ k, b $ are\n",
    "  satisfied.  \n",
    "- $ \\check q = q(K,B) $ and\n",
    "  $ \\check p = p(K,B) $.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4094b58",
   "metadata": {},
   "source": [
    "## Pseudo Code\n",
    "\n",
    "Before displaying our Python code for computing a BCG incomplete markets equilibrium,\n",
    "we’ll sketch some pseudo code that describes its logical flow.\n",
    "\n",
    "Here goes:\n",
    "\n",
    "1. Set upper and lower bounds for firm value as $ V_h $ and\n",
    "  $ V_l $, for capital as $ k_h $ and $ k_l $, and for debt\n",
    "  as $ b_h $ and $ b_l $.  \n",
    "1. Conjecture firm value $ V = \\frac{1}{2}(V_h + V_l) $  \n",
    "1. Conjecture debt level $ b = \\frac{1}{2}(b_h + b_l) $.  \n",
    "1. Conjecture capital $ k = \\frac{1}{2}(k_h + k_l) $.  \n",
    "1. Compute the default threshold\n",
    "  $ \\epsilon^* \\equiv \\log\\left(\\frac{b}{Ak^\\alpha}\\right) $.  \n",
    "1. (In this step we abuse notation by freezing $ V, k, b $ and in\n",
    "  effect temporarily treating them as Big $ K,B $ values. Thus, in\n",
    "  this step 6 little $ k, b $ are frozen at guessed at value of $ K, B $.)\n",
    "  Fixing the values of $ V $, $ b $ and $ k $, compute\n",
    "  optimal choices of consumption $ c^i $ with consumers’ FOCs.\n",
    "  Assume that only agent 2 holds debt: $ \\xi^2 = b $ and that both agents\n",
    "  hold equity: $ 0 <\\theta^i < 1 $ for $ i=1,2 $.  \n",
    "1. Set high and low bounds for equity holdings for agent 1 as  $ \\theta^1_h $ and $ \\theta^1_l $. Guess\n",
    "  $ \\theta^1 = \\frac{1}{2}(\\theta^1_h + \\theta^1_l) $, and\n",
    "  $ \\theta^2 = 1 - \\theta^1 $. While\n",
    "  $ |\\theta^1_h - \\theta^1_l| $ is large:  \n",
    "  - Compute agent 1’s valuation of the equity claim with a\n",
    "    fixed-point iteration:  \n",
    "    $ q_1 = \\beta \\int \\frac{u^\\prime(c^1_1(\\epsilon))}{u^\\prime(c^1_0)} d^e(k,b;\\epsilon) g(\\epsilon) \\ d\\epsilon $  \n",
    "    where  \n",
    "    $ c^1_1(\\epsilon) = w^1_1(\\epsilon) + \\theta^1 d^e(k,b;\\epsilon) $  \n",
    "    and  \n",
    "    $ c^1_0 = w^1_0 + \\theta^1_0V - q_1\\theta^1 $  \n",
    "  - Compute agent 2’s valuation of the bond claim with a\n",
    "    fixed-point iteration:  \n",
    "    $ p = \\beta \\int \\frac{u^\\prime(c^2_1(\\epsilon))}{u^\\prime(c^2_0)} d^b(k,b;\\epsilon) g(\\epsilon) \\ d\\epsilon $  \n",
    "    where  \n",
    "    $ c^2_1(\\epsilon) = w^2_1(\\epsilon) + \\theta^2 d^e(k,b;\\epsilon) + b $  \n",
    "    and  \n",
    "    $ c^2_0 = w^2_0 + \\theta^2_0 V - q_1 \\theta^2 - pb $  \n",
    "  - Compute agent 2’s valuation of the equity claim with a\n",
    "    fixed-point iteration:  \n",
    "    $ q_2 = \\beta \\int \\frac{u^\\prime(c^2_1(\\epsilon))}{u^\\prime(c^2_0)} d^e(k,b;\\epsilon) g(\\epsilon) \\ d\\epsilon $  \n",
    "    where  \n",
    "    $ c^2_1(\\epsilon) = w^2_1(\\epsilon) + \\theta^2 d^e(k,b;\\epsilon) + b $  \n",
    "    and  \n",
    "    $ c^2_0 = w^2_0 + \\theta^2_0 V - q_2 \\theta^2 - pb $  \n",
    "  - If $ q_1 > q_2 $, Set $ \\theta_l = \\theta^1 $;\n",
    "    otherwise, set $ \\theta_h = \\theta^1 $.  \n",
    "  - Repeat steps 6Aa through 6Ad until\n",
    "    $ |\\theta^1_h - \\theta^1_l| $ is small.  \n",
    "1. Set bond price as $ p $ and equity price as  $ q = \\max(q_1,q_2) $.  \n",
    "1. Compute optimal choices of consumption:  \n",
    "  $$\n",
    "  \\begin{aligned}\n",
    "     c^1_0 &= w^1_0 + \\theta^1_0V - q\\theta^1 \\\\\n",
    "     c^2_0 &= w^2_0 + \\theta^2_0V - q\\theta^2 - pb \\\\\n",
    "     c^1_1(\\epsilon) &= w^1_1(\\epsilon) + \\theta^1 d^e(k,b;\\epsilon) \\\\\n",
    "     c^2_1(\\epsilon) &= w^2_1(\\epsilon) + \\theta^2 d^e(k,b;\\epsilon) + b\n",
    "     \\end{aligned}\n",
    "  $$\n",
    "1. (Here we confess to abusing notation again, but now in a different\n",
    "  way. In step 7, we interpret frozen $ c^i $s as Big\n",
    "  $ C^i $. We do this to solve the firm’s problem.) Fixing the\n",
    "  values of $ c^i_0 $ and $ c^i_1(\\epsilon) $, compute optimal\n",
    "  choices of capital $ k $ and debt level $ b $ using the\n",
    "  firm’s first order necessary conditions.  \n",
    "1. Compute deviations from the firm’s FONC for capital $ k $ as:  \n",
    "  $ kfoc = \\beta \\alpha A k^{\\alpha - 1} \\left( \\int \\frac{u^\\prime(c^2_1(\\epsilon))}{u^\\prime(c^2_0)}  e^\\epsilon g(\\epsilon) \\ d\\epsilon \\right) - 1 $  \n",
    "  - If $ kfoc > 0 $, Set $ k_l = k $; otherwise, set\n",
    "    $ k_h = k $.  \n",
    "  - Repeat steps 4 through 7A until $ |k_h-k_l| $ is small.  \n",
    "1. Compute deviations from the firm’s FONC for debt level $ b $  as:  \n",
    "  $ bfoc = \\beta \\left[ \\int_{\\epsilon^*}^\\infty \\left( \\frac{u^\\prime(c^1_1(\\epsilon))}{u^\\prime(c^1_0)} \\right) g(\\epsilon) \\ d\\epsilon -  \\int_{\\epsilon^*}^\\infty \\left( \\frac{u^\\prime(c^2_1(\\epsilon))}{u^\\prime(c^2_0)} \\right)  g(\\epsilon) \\ d\\epsilon \\right] $  \n",
    "  - If $ bfoc > 0 $, Set $ b_h = b $; otherwise, set\n",
    "    $ b_l = b $.  \n",
    "  - Repeat steps 3 through 7B until $ |b_h-b_l| $ is small.  \n",
    "1. Given prices $ q $ and $ p $ from step 6, and the firm\n",
    "  choices of $ k $ and $ b $ from step 7, compute the synthetic\n",
    "  firm value:  \n",
    "  $ V_x = -k + q + pb $  \n",
    "  - If $ V_x > V $, then set $ V_l = V $; otherwise, set\n",
    "    $ V_h = V $.  \n",
    "  - Repeat steps 1 through 8 until $ |V_x - V| $ is small.  \n",
    "1. Ultimately, the algorithm returns  equilibrium capital\n",
    "  $ k^* $, debt $ b^* $ and firm value $ V^* $, as well as\n",
    "  the following equilibrium values:  \n",
    "  - Equity holdings $ \\theta^{1,*} = \\theta^1(k^*,b^*) $  \n",
    "  - Prices $ q^*=q(k^*,b^*), \\ p^*=p(k^*,b^*) $  \n",
    "  - Consumption plans\n",
    "    $ C^{1,*}_0 = c^1_0(k^*,b^*),\\ C^{2,*}_0 = c^2_0(k^*,b^*),  \\ C^{1,*}_1(\\epsilon) = c^1_1(k^*,b^*;\\epsilon),\\ C^{1,*}_1(\\epsilon) = c^2_1(k^*,b^*;\\epsilon) $.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fca85145",
   "metadata": {},
   "source": [
    "## Code\n",
    "\n",
    "We create a Python class `BCG_incomplete_markets` to compute the\n",
    "equilibrium allocations of the incomplete market BCG model, given a set\n",
    "of parameter values.\n",
    "\n",
    "The class includes the following methods,  i.e., functions:\n",
    "\n",
    "- `solve_eq`: solves the BCG model and returns the equilibrium values\n",
    "  of capital $ k $, debt $ b $ and firm value $ V $, as\n",
    "  well as  \n",
    "  - agent 1’s equity holdings $ \\theta^{1,*} $  \n",
    "  - prices $ q^*, p^* $  \n",
    "  - consumption plans\n",
    "    $ C^{1,*}_0, C^{2,*}_0, C^{1,*}_1(\\epsilon), C^{2,*}_1(\\epsilon) $.  \n",
    "- `eq_valuation`: inputs equilibrium consumpion plans $ C^* $ and\n",
    "  outputs the following valuations for each pair of $ (k,b) $ in\n",
    "  the grid:  \n",
    "  - the firm $ V(k,b) $  \n",
    "  - the equity $ q(k,b) $  \n",
    "  - the bond $ p(k,b) $.  \n",
    "\n",
    "\n",
    "Parameters include:\n",
    "\n",
    "- $ \\chi_1 $, $ \\chi_2 $:  correlation parameter for agent 1\n",
    "  and 2. Default values are respectively 0 and 0.9.  \n",
    "- $ w^1_0 $, $ w^2_0 $:  initial endowments. Default values\n",
    "  are respectively 0.9 and 1.1.  \n",
    "- $ \\theta^1_0 $, $ \\theta^2_0 $:  initial holding of the\n",
    "  firm. Default values are 0.5.  \n",
    "- $ \\psi $:  risk parameter. Default value is 3.  \n",
    "- $ \\alpha $: Production function parameter. Default value\n",
    "  is 0.6.  \n",
    "- $ A $: Productivity of the firm. Default value is 2.5.  \n",
    "- $ \\mu $, $ \\sigma $: Mean and standard deviation of the\n",
    "  shock distribution. Default values are respectively -0.025 and 0.4  \n",
    "- $ \\beta $: Discount factor. Default value is 0.96.  \n",
    "- bound: Bound for truncated normal distribution. Default value is 3.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0bf31915",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import truncnorm\n",
    "from scipy.integrate import quad\n",
    "from numba import njit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c335a592",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "class BCG_incomplete_markets:\n",
    "\n",
    "    # init method or constructor\n",
    "    def __init__(self,\n",
    "                 𝜒1 = 0,\n",
    "                 𝜒2 = 0.9,\n",
    "                 w10 = 0.9,\n",
    "                 w20 = 1.1,\n",
    "                 𝜃10 = 0.5,\n",
    "                 𝜃20 = 0.5,\n",
    "                 𝜓1 = 3,\n",
    "                 𝜓2 = 3,\n",
    "                 𝛼 = 0.6,\n",
    "                 A = 2.5,\n",
    "                 𝜇 = -0.025,\n",
    "                 𝜎 = 0.4,\n",
    "                 𝛽 = 0.96,\n",
    "                 bound = 3,\n",
    "                 Vl = 0,\n",
    "                 Vh = 0.5,\n",
    "                 kbot = 0.01,\n",
    "                 #ktop = (𝛼*A)**(1/(1-𝛼)),\n",
    "                 ktop = 0.25,\n",
    "                 bbot = 0.1,\n",
    "                 btop = 0.8):\n",
    "\n",
    "        #=========== Setup ===========#\n",
    "        # Risk parameters\n",
    "        self.𝜒1 = 𝜒1\n",
    "        self.𝜒2 = 𝜒2\n",
    "\n",
    "        # Other parameters\n",
    "        self.𝜓1 = 𝜓1\n",
    "        self.𝜓2 = 𝜓2\n",
    "        self.𝛼 = 𝛼\n",
    "        self.A = A\n",
    "        self.𝜇 = 𝜇\n",
    "        self.𝜎 = 𝜎\n",
    "        self.𝛽 = 𝛽\n",
    "        self.bound = bound\n",
    "\n",
    "        # Bounds for firm value, capital, and debt\n",
    "        self.Vl = Vl\n",
    "        self.Vh = Vh\n",
    "        self.kbot = kbot\n",
    "        #self.kbot = (𝛼*A)**(1/(1-𝛼))\n",
    "        self.ktop = ktop\n",
    "        self.bbot = bbot\n",
    "        self.btop = btop\n",
    "\n",
    "        # Utility\n",
    "        self.u = njit(lambda c: (c**(1-𝜓)) / (1-𝜓))\n",
    "\n",
    "        # Initial endowments\n",
    "        self.w10 = w10\n",
    "        self.w20 = w20\n",
    "        self.w0 = w10 + w20\n",
    "\n",
    "        # Initial holdings\n",
    "        self.𝜃10 = 𝜃10\n",
    "        self.𝜃20 = 𝜃20\n",
    "\n",
    "        # Endowments at t=1\n",
    "        self.w11 = njit(lambda 𝜖: np.exp(-𝜒1*𝜇 - 0.5*(𝜒1**2)*(𝜎**2) + 𝜒1*𝜖))\n",
    "        self.w21 = njit(lambda 𝜖: np.exp(-𝜒2*𝜇 - 0.5*(𝜒2**2)*(𝜎**2) + 𝜒2*𝜖))\n",
    "        self.w1 = njit(lambda 𝜖: self.w11(𝜖) + self.w21(𝜖))\n",
    "\n",
    "        # Truncated normal\n",
    "        ta, tb = (-bound - 𝜇) / 𝜎, (bound - 𝜇) / 𝜎\n",
    "        rv = truncnorm(ta, tb, loc=𝜇, scale=𝜎)\n",
    "        𝜖_range = np.linspace(ta, tb, 1000000)\n",
    "        pdf_range = rv.pdf(𝜖_range)\n",
    "        self.g = njit(lambda 𝜖: np.interp(𝜖, 𝜖_range, pdf_range))\n",
    "\n",
    "\n",
    "    #*************************************************************\n",
    "    # Function: Solve for equilibrium of the BCG model\n",
    "    #*************************************************************\n",
    "    def solve_eq(self, print_crit=True):\n",
    "\n",
    "        # Load parameters\n",
    "        𝜓1 = self.𝜓1\n",
    "        𝜓2 = self.𝜓2\n",
    "        𝛼 = self.𝛼\n",
    "        A = self.A\n",
    "        𝛽 = self.𝛽\n",
    "        bound = self.bound\n",
    "        Vl = self.Vl\n",
    "        Vh = self.Vh\n",
    "        kbot = self.kbot\n",
    "        ktop = self.ktop\n",
    "        bbot = self.bbot\n",
    "        btop = self.btop\n",
    "        w10 = self.w10\n",
    "        w20 = self.w20\n",
    "        𝜃10 = self.𝜃10\n",
    "        𝜃20 = self.𝜃20\n",
    "        w11 = self.w11\n",
    "        w21 = self.w21\n",
    "        g = self.g\n",
    "\n",
    "        # We need to find a fixed point on the value of the firm\n",
    "        V_crit = 1\n",
    "\n",
    "        Y = njit(lambda 𝜖, fk: np.exp(𝜖)*fk)\n",
    "        intqq1 = njit(lambda 𝜖, fk, 𝜃1, 𝜓1, b: (w11(𝜖) + 𝜃1*(Y(𝜖, fk) - b))**(-𝜓1)*(Y(𝜖, fk) - b)*g(𝜖))\n",
    "        intp1 = njit(lambda 𝜖, fk, 𝜓2, b: (Y(𝜖, fk)/b)*(w21(𝜖) + Y(𝜖, fk))**(-𝜓2)*g(𝜖))\n",
    "        intp2 = njit(lambda 𝜖, fk, 𝜃2, 𝜓2, b: (w21(𝜖) + 𝜃2*(Y(𝜖, fk)-b) + b)**(-𝜓2)*g(𝜖))\n",
    "        intqq2 = njit(lambda 𝜖, fk, 𝜃2, 𝜓2, b: (w21(𝜖) + 𝜃2*(Y(𝜖, fk)-b) + b)**(-𝜓2)*(Y(𝜖, fk) - b)*g(𝜖))\n",
    "        intk1 = njit(lambda 𝜖, fk, 𝜓2: (w21(𝜖) + Y(𝜖, fk))**(-𝜓2)*np.exp(𝜖)*g(𝜖))\n",
    "        intk2 = njit(lambda 𝜖, fk, 𝜃2, 𝜓2, b: (w21(𝜖) + 𝜃2*(Y(𝜖, fk)-b) + b)**(-𝜓2)*np.exp(𝜖)*g(𝜖))\n",
    "        intB1 = njit(lambda 𝜖, fk, 𝜃1, 𝜓1, b: (w11(𝜖) + 𝜃1*(Y(𝜖, fk) - b))**(-𝜓1)*g(𝜖))\n",
    "        intB2 = njit(lambda 𝜖, fk, 𝜃2, 𝜓2, b: (w21(𝜖) + 𝜃2*(Y(𝜖, fk) - b) + b)**(-𝜓2)*g(𝜖))\n",
    "\n",
    "        while V_crit>1e-4:\n",
    "\n",
    "            # We begin by adding the guess for the value of the firm to endowment\n",
    "            V = (Vl+Vh)/2\n",
    "            ww10 = w10 + 𝜃10*V\n",
    "            ww20 = w20 + 𝜃20*V\n",
    "\n",
    "            # Figure out the optimal level of debt\n",
    "            bl = bbot\n",
    "            bh = btop\n",
    "            b_crit=1\n",
    "\n",
    "            while b_crit>1e-5:\n",
    "\n",
    "                # Setting the conjecture for debt\n",
    "                b = (bl+bh)/2\n",
    "\n",
    "                # Figure out the optimal level of capital\n",
    "                kl = kbot\n",
    "                kh = ktop\n",
    "                k_crit=1\n",
    "\n",
    "                while k_crit>1e-5:\n",
    "\n",
    "                    # Setting the conjecture for capital\n",
    "                    k = (kl+kh)/2\n",
    "\n",
    "                    # Production\n",
    "                    fk = A*(k**𝛼)\n",
    "#                     Y = lambda 𝜖: np.exp(𝜖)*fk\n",
    "\n",
    "                    # Compute integration threshold\n",
    "                    epstar = np.log(b/fk)\n",
    "\n",
    "\n",
    "                    #**************************************************************\n",
    "                    # Compute the prices and allocations consistent with consumers'\n",
    "                    # Euler equations\n",
    "                    #**************************************************************\n",
    "\n",
    "                    # We impose the following:\n",
    "                    # Agent 1 buys equity\n",
    "                    # Agent 2 buys equity and all debt\n",
    "                    # Agents trade such that prices converge\n",
    "\n",
    "                    #========\n",
    "                    # Agent 1\n",
    "                    #========\n",
    "                    # Holdings\n",
    "                    𝜉1 = 0\n",
    "                    𝜃1a = 0.3\n",
    "                    𝜃1b = 1\n",
    "\n",
    "                    while abs(𝜃1b - 𝜃1a) > 0.001:\n",
    "\n",
    "                        𝜃1 = (𝜃1a + 𝜃1b) / 2\n",
    "\n",
    "                        # qq1 is the equity price consistent with agent-1 Euler Equation\n",
    "                        ## Note: Price is in the date-0 budget constraint of the agent\n",
    "\n",
    "                        ## First, compute the constant term that is not influenced by q\n",
    "                        ## that is, 𝛽E[u'(c^{1}_{1})d^{e}(k,B)]\n",
    "#                         intqq1 = lambda 𝜖: (w11(𝜖) + 𝜃1*(Y(𝜖, fk) - b))**(-𝜓1)*(Y(𝜖, fk) - b)*g(𝜖)\n",
    "#                         const_qq1 = 𝛽 * quad(intqq1,epstar,bound)[0]\n",
    "\n",
    "                        const_qq1 = 𝛽 * quad(intqq1,epstar,bound, args=(fk, 𝜃1, 𝜓1, b))[0]\n",
    "\n",
    "\n",
    "                        ## Second, iterate to get the equity price q\n",
    "                        qq1l = 0\n",
    "                        qq1h = ww10\n",
    "                        diff = 1\n",
    "                        while diff > 1e-7:\n",
    "                            qq1 = (qq1l+qq1h)/2\n",
    "                            rhs = const_qq1/((ww10-qq1*𝜃1)**(-𝜓1));\n",
    "                            if (rhs > qq1):\n",
    "                                qq1l = qq1\n",
    "                            else:\n",
    "                                qq1h = qq1\n",
    "                            diff = abs(qq1l-qq1h)\n",
    "\n",
    "                        #========\n",
    "                        # Agent 2\n",
    "                        #========\n",
    "                        𝜉2 = b - 𝜉1\n",
    "                        𝜃2 = 1 - 𝜃1\n",
    "\n",
    "                        # p is the bond price consistent with agent-2 Euler Equation\n",
    "                        ## Note: Price is in the date-0 budget constraint of the agent\n",
    "\n",
    "                        ## First, compute the constant term that is not influenced by p\n",
    "                        ## that is, 𝛽E[u'(c^{2}_{1})d^{b}(k,B)]\n",
    "#                         intp1 = lambda 𝜖: (Y(𝜖, fk)/b)*(w21(𝜖) + Y(𝜖, fk))**(-𝜓2)*g(𝜖)\n",
    "#                         intp2 = lambda 𝜖: (w21(𝜖) + 𝜃2*(Y(𝜖, fk)-b) + b)**(-𝜓2)*g(𝜖)\n",
    "#                         const_p = 𝛽 * (quad(intp1,-bound,epstar)[0] + quad(intp2,epstar,bound)[0])\n",
    "                        const_p = 𝛽 * (quad(intp1,-bound,epstar, args=(fk, 𝜓2, b))[0]\\\n",
    "                                        + quad(intp2,epstar,bound, args=(fk, 𝜃2, 𝜓2, b))[0])\n",
    "\n",
    "                        ## iterate to get the bond price p\n",
    "                        pl = 0\n",
    "                        ph = ww20/b\n",
    "                        diff = 1\n",
    "                        while diff > 1e-7:\n",
    "                            p = (pl+ph)/2\n",
    "                            rhs = const_p/((ww20-qq1*𝜃2-p*b)**(-𝜓2))\n",
    "                            if (rhs > p):\n",
    "                                pl = p\n",
    "                            else:\n",
    "                                ph = p\n",
    "                            diff = abs(pl-ph)\n",
    "\n",
    "                        # qq2 is the equity price consistent with agent-2 Euler Equation\n",
    "#                         intqq2 = lambda 𝜖: (w21(𝜖) + 𝜃2*(Y(𝜖, fk)-b) + b)**(-𝜓2)*(Y(𝜖, fk) - b)*g(𝜖)\n",
    "                        const_qq2 = 𝛽 * quad(intqq2,epstar,bound, args=(fk, 𝜃2, 𝜓2, b))[0]\n",
    "                        qq2l = 0\n",
    "                        qq2h = ww20\n",
    "                        diff = 1\n",
    "                        while diff > 1e-7:\n",
    "                            qq2 = (qq2l+qq2h)/2\n",
    "                            rhs = const_qq2/((ww20-qq2*𝜃2-p*b)**(-𝜓2));\n",
    "                            if (rhs > qq2):\n",
    "                                qq2l = qq2\n",
    "                            else:\n",
    "                                qq2h = qq2\n",
    "                            diff = abs(qq2l-qq2h)\n",
    "\n",
    "                        # q be the maximum valuation for the equity among agents\n",
    "                        ## This will be the equity price based on Makowski's criterion\n",
    "                        q = max(qq1,qq2)\n",
    "\n",
    "                        #================\n",
    "                        # Update holdings\n",
    "                        #================\n",
    "                        if qq1 > qq2:\n",
    "                            𝜃1a = 𝜃1\n",
    "                        else:\n",
    "                            𝜃1b = 𝜃1\n",
    "\n",
    "                    #================\n",
    "                    # Get consumption\n",
    "                    #================\n",
    "                    c10 = ww10 - q*𝜃1\n",
    "                    c11 = lambda 𝜖: w11(𝜖) + 𝜃1*max(Y(𝜖, fk)-b,0)\n",
    "                    c20 = ww20 - q*(1-𝜃1) - p*b\n",
    "                    c21 = lambda 𝜖: w21(𝜖) + (1-𝜃1)*max(Y(𝜖, fk)-b,0) + min(Y(𝜖, fk),b)\n",
    "\n",
    "\n",
    "                    #*************************************************\n",
    "                    # Compute the first order conditions for the firm\n",
    "                    #*************************************************\n",
    "\n",
    "                    #===========\n",
    "                    # Equity FOC\n",
    "                    #===========\n",
    "                    # Only agent 2's IMRS is relevent\n",
    "#                     intk1 = lambda 𝜖: (w21(𝜖) + Y(𝜖, fk))**(-𝜓2)*np.exp(𝜖)*g(𝜖)\n",
    "#                     intk2 = lambda 𝜖: (w21(𝜖) + 𝜃2*(Y(𝜖, fk)-b) + b)**(-𝜓2)*np.exp(𝜖)*g(𝜖)\n",
    "#                     kfoc_num = quad(intk1,-bound,epstar)[0] + quad(intk2,epstar,bound)[0]\n",
    "                    kfoc_num = quad(intk1,-bound,epstar, args=(fk, 𝜓2))[0] + quad(intk2,epstar,bound, args=(fk, 𝜃2, 𝜓2, b))[0]\n",
    "                    kfoc_denom = (ww20- q*𝜃2 - p*b)**(-𝜓2)\n",
    "                    kfoc = 𝛽*𝛼*A*(k**(𝛼-1))*(kfoc_num/kfoc_denom) - 1\n",
    "\n",
    "                    if (kfoc > 0):\n",
    "                        kl = k\n",
    "                    else:\n",
    "                        kh = k\n",
    "                    k_crit = abs(kh-kl)\n",
    "\n",
    "                    if print_crit:\n",
    "                        print(\"critical value of k: {:.5f}\".format(k_crit))\n",
    "\n",
    "\n",
    "                #=========\n",
    "                # Bond FOC\n",
    "                #=========\n",
    "#                 intB1 = lambda 𝜖: (w11(𝜖) + 𝜃1*(Y(𝜖, fk) - b))**(-𝜓1)*g(𝜖)\n",
    "#                 intB2 = lambda 𝜖: (w21(𝜖) + 𝜃2*(Y(𝜖, fk) - b) + b)**(-𝜓2)*g(𝜖)\n",
    "\n",
    "#                 bfoc1 = quad(intB1,epstar,bound)[0] / (ww10 - q*𝜃1)**(-𝜓1)\n",
    "#                 bfoc2 = quad(intB2,epstar,bound)[0] / (ww20 - q*𝜃2 - p*b)**(-𝜓2)\n",
    "\n",
    "                bfoc1 = quad(intB1,epstar,bound, args=(fk, 𝜃1, 𝜓1, b))[0] / (ww10 - q*𝜃1)**(-𝜓1)\n",
    "                bfoc2 = quad(intB2,epstar,bound, args=(fk, 𝜃2, 𝜓2, b))[0] / (ww20 - q*𝜃2 - p*b)**(-𝜓2)\n",
    "                bfoc = bfoc1 - bfoc2\n",
    "\n",
    "                if (bfoc > 0):\n",
    "                    bh = b\n",
    "                else:\n",
    "                    bl = b\n",
    "                b_crit = abs(bh-bl)\n",
    "\n",
    "                if print_crit:\n",
    "                    print(\"#=== critical value of b: {:.5f}\".format(b_crit))\n",
    "\n",
    "            # Compute the value of the firm\n",
    "            value_x = -k + q + p*b\n",
    "            if (value_x > V):\n",
    "                 Vl = V\n",
    "            else:\n",
    "                 Vh = V\n",
    "            V_crit = abs(value_x-V)\n",
    "\n",
    "            if print_crit:\n",
    "                print(\"#====== critical value of V: {:.5f}\".format(V_crit))\n",
    "\n",
    "            print('k,b,p,q,kfoc,bfoc,epstar,V,V_crit')\n",
    "            formattedList = [\"%.3f\" % member for member in [k,\n",
    "                                                            b,\n",
    "                                                            p,\n",
    "                                                            q,\n",
    "                                                            kfoc,\n",
    "                                                            bfoc,\n",
    "                                                            epstar,\n",
    "                                                            V,\n",
    "                                                            V_crit]]\n",
    "            print(formattedList)\n",
    "\n",
    "        #*********************************\n",
    "        # Equilibrium values\n",
    "        #*********************************\n",
    "\n",
    "        # Return the results\n",
    "        kss = k\n",
    "        bss = b\n",
    "        Vss = V\n",
    "        qss = q\n",
    "        pss = p\n",
    "        c10ss = c10\n",
    "        c11ss = c11\n",
    "        c20ss = c20\n",
    "        c21ss = c21\n",
    "        𝜃1ss = 𝜃1\n",
    "\n",
    "\n",
    "        # Print the results\n",
    "        print('finished')\n",
    "       # print('k,b,p,q,kfoc,bfoc,epstar,V,V_crit')\n",
    "        #formattedList = [\"%.3f\" % member for member in [kss,\n",
    "        #                                                bss,\n",
    "        #                                                pss,\n",
    "        #                                                qss,\n",
    "        #                                                kfoc,\n",
    "        #                                                bfoc,\n",
    "        #                                                epstar,\n",
    "        #                                                Vss,\n",
    "        #                                                V_crit]]\n",
    "        #print(formattedList)\n",
    "\n",
    "        return kss,bss,Vss,qss,pss,c10ss,c11ss,c20ss,c21ss,𝜃1ss\n",
    "\n",
    "\n",
    "    #*************************************************************\n",
    "    # Function: Equity and bond valuations by different agents\n",
    "    #*************************************************************\n",
    "    def valuations_by_agent(self,\n",
    "                            c10, c11, c20, c21,\n",
    "                            k, b):\n",
    "\n",
    "        # Load parameters\n",
    "        𝜓1 = self.𝜓1\n",
    "        𝜓2 = self.𝜓2\n",
    "        𝛼 = self.𝛼\n",
    "        A = self.A\n",
    "        𝛽 = self.𝛽\n",
    "        bound = self.bound\n",
    "        Vl = self.Vl\n",
    "        Vh = self.Vh\n",
    "        kbot = self.kbot\n",
    "        ktop = self.ktop\n",
    "        bbot = self.bbot\n",
    "        btop = self.btop\n",
    "        w10 = self.w10\n",
    "        w20 = self.w20\n",
    "        𝜃10 = self.𝜃10\n",
    "        𝜃20 = self.𝜃20\n",
    "        w11 = self.w11\n",
    "        w21 = self.w21\n",
    "        g = self.g\n",
    "\n",
    "        # Get functions for IMRS/state price density\n",
    "        IMRS1 = lambda 𝜖: 𝛽 * (c11(𝜖)/c10)**(-𝜓1)*g(𝜖)\n",
    "        IMRS2 = lambda 𝜖: 𝛽 * (c21(𝜖)/c20)**(-𝜓2)*g(𝜖)\n",
    "\n",
    "        # Production\n",
    "        fk = A*(k**𝛼)\n",
    "        Y = lambda 𝜖: np.exp(𝜖)*fk\n",
    "\n",
    "        # Compute integration threshold\n",
    "        epstar = np.log(b/fk)\n",
    "\n",
    "        # Compute equity valuation with agent 1's IMRS\n",
    "        intQ1 = lambda 𝜖: IMRS1(𝜖)*(Y(𝜖) - b)\n",
    "        Q1 = quad(intQ1, epstar, bound)[0]\n",
    "\n",
    "        # Compute bond valuation with agent 1's IMRS\n",
    "        intP1 = lambda 𝜖: IMRS1(𝜖)*Y(𝜖)/b\n",
    "        P1 = quad(intP1, -bound, epstar)[0] + quad(IMRS1, epstar, bound)[0]\n",
    "\n",
    "        # Compute equity valuation with agent 2's IMRS\n",
    "        intQ2 = lambda 𝜖: IMRS2(𝜖)*(Y(𝜖) - b)\n",
    "        Q2 = quad(intQ2, epstar, bound)[0]\n",
    "\n",
    "        # Compute bond valuation with agent 2's IMRS\n",
    "        intP2 = lambda 𝜖: IMRS2(𝜖)*Y(𝜖)/b\n",
    "        P2 = quad(intP2, -bound, epstar)[0] + quad(IMRS2, epstar, bound)[0]\n",
    "\n",
    "        return Q1,Q2,P1,P2\n",
    "\n",
    "\n",
    "    #*************************************************************\n",
    "    # Function: equilibrium valuations for firm, equity, bond\n",
    "    #*************************************************************\n",
    "    def eq_valuation(self, c10, c11, c20, c21, N=30):\n",
    "\n",
    "        # Load parameters\n",
    "        𝜓1 = self.𝜓1\n",
    "        𝜓2 = self.𝜓2\n",
    "        𝛼 = self.𝛼\n",
    "        A = self.A\n",
    "        𝛽 = self.𝛽\n",
    "        bound = self.bound\n",
    "        Vl = self.Vl\n",
    "        Vh = self.Vh\n",
    "        kbot = self.kbot\n",
    "        ktop = self.ktop\n",
    "        bbot = self.bbot\n",
    "        btop = self.btop\n",
    "        w10 = self.w10\n",
    "        w20 = self.w20\n",
    "        𝜃10 = self.𝜃10\n",
    "        𝜃20 = self.𝜃20\n",
    "        w11 = self.w11\n",
    "        w21 = self.w21\n",
    "        g = self.g\n",
    "\n",
    "        # Create grids\n",
    "        kgrid, bgrid = np.meshgrid(np.linspace(kbot,ktop,N),\n",
    "                                   np.linspace(bbot,btop,N))\n",
    "        Vgrid = np.zeros_like(kgrid)\n",
    "        Qgrid = np.zeros_like(kgrid)\n",
    "        Pgrid = np.zeros_like(kgrid)\n",
    "\n",
    "        # Loop: firm value\n",
    "        for i in range(N):\n",
    "            for j in range(N):\n",
    "\n",
    "                # Get capital and debt\n",
    "                k = kgrid[i,j]\n",
    "                b = bgrid[i,j]\n",
    "\n",
    "                # Valuations by each agent\n",
    "                Q1,Q2,P1,P2 = self.valuations_by_agent(c10,\n",
    "                                                       c11,\n",
    "                                                       c20,\n",
    "                                                       c21,\n",
    "                                                       k,\n",
    "                                                       b)\n",
    "\n",
    "                # The prices will be the maximum of the valuations\n",
    "                Q = max(Q1,Q2)\n",
    "                P = max(P1,P2)\n",
    "\n",
    "                # Compute firm value\n",
    "                V = -k + Q + P*b\n",
    "                Vgrid[i,j] = V\n",
    "                Qgrid[i,j] = Q\n",
    "                Pgrid[i,j] = P\n",
    "\n",
    "        return kgrid, bgrid, Vgrid, Qgrid, Pgrid"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2be2452e",
   "metadata": {},
   "source": [
    "## Examples\n",
    "\n",
    "Below we show some examples computed with the class `BCG_incomplete markets`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46618628",
   "metadata": {},
   "source": [
    "### First example\n",
    "\n",
    "In the first example, we set up an instance of the BCG incomplete\n",
    "markets model with default parameter values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5869bbb5",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "mdl = BCG_incomplete_markets()\n",
    "kss,bss,Vss,qss,pss,c10ss,c11ss,c20ss,c21ss,𝜃1ss = mdl.solve_eq(print_crit=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a2a26ad",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "print(-kss+qss+pss*bss)\n",
    "print(Vss)\n",
    "print(𝜃1ss)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a14b5cb9",
   "metadata": {},
   "source": [
    "Python reports to us that the equilibrium firm value is $ V=0.101 $,\n",
    "with capital $ k = 0.151 $ and debt $ b=0.484 $.\n",
    "\n",
    "Let’s verify some things that have to be true if our algorithm has truly\n",
    "found an equilibrium.\n",
    "\n",
    "Thus, let’s see if the firm is actually maximizing its firm value given\n",
    "the equilibrium pricing function $ q(k,b) $ for equity and\n",
    "$ p(k,b) $ for  bonds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb7e7489",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "kgrid, bgrid, Vgrid, Qgrid, Pgrid = mdl.eq_valuation(c10ss, c11ss, c20ss, c21ss,N=30)\n",
    "\n",
    "print('Maximum valuation of the firm value in the (k,B) grid: {:.5f}'.format(Vgrid.max()))\n",
    "print('Equilibrium firm value: {:.5f}'.format(Vss))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f77c836",
   "metadata": {},
   "source": [
    "Up to the approximation involved in using a discrete grid, these numbers\n",
    "give us comfort that the firm does indeed seem to be maximizing its\n",
    "value at the top of the value hill on the $ (k,b) $ plane that it\n",
    "faces.\n",
    "\n",
    "Below we will plot the firm’s value as a function of $ k,b $.\n",
    "\n",
    "We’ll also plot the equilibrium price functions $ q(k,b) $ and\n",
    "$ p(k,b) $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f3a31b8",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from IPython.display import Image\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits import mplot3d\n",
    "import plotly.graph_objs as go\n",
    "\n",
    "# Firm Valuation\n",
    "fig = go.Figure(data=[go.Scatter3d(x=[kss],\n",
    "                                   y=[bss],\n",
    "                                   z=[Vss],\n",
    "                                   mode='markers',\n",
    "                                   marker=dict(size=3, color='red')),\n",
    "                      go.Surface(x=kgrid,\n",
    "                                 y=bgrid,\n",
    "                                 z=Vgrid,\n",
    "                                 colorscale='Greens',opacity=0.6)])\n",
    "\n",
    "fig.update_layout(scene = dict(\n",
    "                    xaxis_title='x - Capital k',\n",
    "                    yaxis_title='y - Debt b',\n",
    "                    zaxis_title='z - Firm Value V',\n",
    "                    aspectratio = dict(x=1,y=1,z=1)),\n",
    "                  width=700,\n",
    "                  height=700,\n",
    "                  margin=dict(l=50, r=50, b=65, t=90))\n",
    "fig.update_layout(scene_camera=dict(eye=dict(x=1.5, y=-1.5, z=2)))\n",
    "fig.update_layout(title='Equilibrium firm valuation for the grid of (k,b)')\n",
    "\n",
    "# Export to PNG file\n",
    "Image(fig.to_image(format=\"png\"))\n",
    "# fig.show() will provide interactive plot when running\n",
    "# code locally"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7696c0fa",
   "metadata": {},
   "source": [
    "#### A Modigliani-Miller theorem?\n",
    "\n",
    "The red dot in the above graph is **both** an equilibrium $ (b,k) $\n",
    "chosen by a representative firm **and** the equilibrium $ B, K $\n",
    "pair chosen by the aggregate of all firms.\n",
    "\n",
    "Thus, **in equilibrium** it\n",
    "is true that\n",
    "\n",
    "$$\n",
    "(b,k) = (B,K)\n",
    "$$\n",
    "\n",
    "But an individual firm named $ \\zeta \\in [0,1] $ neither knows nor\n",
    "cares whether it sets $ (b(\\zeta),k(\\zeta)) = (B,K) $.\n",
    "\n",
    "Indeed the above graph has a ridge of $ b(\\zeta) $’s that also\n",
    "maximize the firm’s value so long as it sets $ k(\\zeta) = K $.\n",
    "\n",
    "Here it is important that the measure of firms that deviate from setting\n",
    "$ b $ at the red dot is very small – measure zero – so that\n",
    "$ B $ remains at the red dot even while one firm $ \\zeta $\n",
    "deviates.\n",
    "\n",
    "So within this equilibrium, there is a  *qualified* Modigliani-Miller theorem\n",
    "that asserts that firm $ \\zeta $’s value is\n",
    "independent of how it mixes its financing between equity and bonds (so\n",
    "long as it is not  what other firms do on average).\n",
    "\n",
    "Thus, while an individual firm $ \\zeta $’s financial structure is\n",
    "indeterminate, the **market’s** financial structure is determinant and\n",
    "sits at the red dot in the above graph.\n",
    "\n",
    "This contrasts sharply with the *unqualified* Modigliani-Miller theorem\n",
    "descibed in the complete markets model in the lecture [Irrelevance of Capital Structure with Complete Markets](https://python-advanced.quantecon.org/BCG_complete_mkts.html).\n",
    "\n",
    "There the **market’s** financial structure was indeterminate.\n",
    "\n",
    "These subtle distinctions  bear more thought and exploration.\n",
    "\n",
    "So we will do some calculations  to ferret out a sense in which\n",
    "the equilibrium $ (k,b) = (K,B) $ outcome at the red dot in the\n",
    "above graph is **stable**.\n",
    "\n",
    "In particular, we’ll explore the consequences of some choices of\n",
    "$ b=B $ that deviate from the red dot and ask whether firm\n",
    "$ \\zeta $ would want to remain at that $ b $.\n",
    "\n",
    "In more detail, here is what we’ll do:\n",
    "\n",
    "1. Obtain equilibrium values of capital and debt as $ k^*=K $ and\n",
    "  $ b^*=B $,  the red dot above.  \n",
    "1. Now fix $ k^* $ and let $ b^{**} = b^* - e $ for some\n",
    "  $ e > 0 $. Conjecture that big $ K = k^* $ but big\n",
    "  $ B = b^{**} $.  \n",
    "1. Take $ K $ and $ B $ and compute intertermporal marginal rates of substitution (IMRS’s) as we did before.  \n",
    "1. Taking the **new** IMRS to the firm’s problem. Plot 3D surface for\n",
    "  the valuations of the firm with this **new** IMRS.  \n",
    "1. Check if the value at $ k^* $, $ b^{**} $ is at the top of\n",
    "  this new 3D surface.  \n",
    "1. Repeat these calculations for $ b^{**} = b^* + e $.  \n",
    "\n",
    "\n",
    "To conduct the above procedures, we create a function `off_eq_check`\n",
    "that inputs the BCG model instance parameters, equilibrium capital\n",
    "$ K=k^* $ and debt $ B=b^* $, and a perturbation of debt $ e $.\n",
    "\n",
    "The function outputs the fixed point firm values $ V^{**} $, prices\n",
    "$ q^{**} $, $ p^{**} $, and consumption choices $ c^{**} $.\n",
    "\n",
    "Importantly, we relax the condition that only agent 2 holds bonds.\n",
    "\n",
    "Now **both** agents can hold bonds, i.e., $ 0\\leq \\xi^1 \\leq B $ and\n",
    "$ \\xi^1 +\\xi^2 = B $.\n",
    "\n",
    "That implies the consumers’ budget constraints are:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "c^1_0 &= w^1_0 + \\theta^1_0V - q\\theta^1 - p\\xi^1 \\\\\n",
    "c^2_0 &= w^2_0 + \\theta^2_0V - q\\theta^2 - p\\xi^2 \\\\\n",
    "c^1_1(\\epsilon) &= w^1_1(\\epsilon) + \\theta^1 d^e(k,b;\\epsilon) + \\xi^1 \\\\\n",
    "c^2_1(\\epsilon) &= w^2_1(\\epsilon) + \\theta^2 d^e(k,b;\\epsilon) + \\xi^2\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The function also outputs agent 1’s bond holdings $ \\xi_1 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f45a5f26",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def off_eq_check(mdl,kss,bss,e=0.1):\n",
    "    # Big K and big B\n",
    "    k = kss\n",
    "    b = bss + e\n",
    "\n",
    "    # Load parameters\n",
    "    𝜓1 = mdl.𝜓1\n",
    "    𝜓2 = mdl.𝜓2\n",
    "    𝛼 = mdl.𝛼\n",
    "    A = mdl.A\n",
    "    𝛽 = mdl.𝛽\n",
    "    bound = mdl.bound\n",
    "    Vl = mdl.Vl\n",
    "    Vh = mdl.Vh\n",
    "    kbot = mdl.kbot\n",
    "    ktop = mdl.ktop\n",
    "    bbot = mdl.bbot\n",
    "    btop = mdl.btop\n",
    "    w10 = mdl.w10\n",
    "    w20 = mdl.w20\n",
    "    𝜃10 = mdl.𝜃10\n",
    "    𝜃20 = mdl.𝜃20\n",
    "    w11 = mdl.w11\n",
    "    w21 = mdl.w21\n",
    "    g = mdl.g\n",
    "\n",
    "    Y = njit(lambda 𝜖, fk: np.exp(𝜖)*fk)\n",
    "    intqq1 = njit(lambda 𝜖, fk, 𝜃1, 𝜓1, 𝜉1, b: (w11(𝜖) + 𝜃1*(Y(𝜖, fk) - b) + 𝜉1)**(-𝜓1)*(Y(𝜖, fk) - b)*g(𝜖))\n",
    "    intpp1a = njit(lambda 𝜖, fk, 𝜓1, 𝜉1, b: (Y(𝜖, fk)/b)*(w11(𝜖) + Y(𝜖, fk)/b*𝜉1)**(-𝜓1)*g(𝜖))\n",
    "    intpp1b = njit(lambda 𝜖, fk, 𝜃1, 𝜓1, 𝜉1, b: (w11(𝜖) + 𝜃1*(Y(𝜖, fk)-b) + 𝜉1)**(-𝜓1)*g(𝜖))\n",
    "    intpp2a = njit(lambda 𝜖, fk, 𝜓2, 𝜉2, b: (Y(𝜖, fk)/b)*(w21(𝜖) + Y(𝜖, fk)/b*𝜉2)**(-𝜓2)*g(𝜖))\n",
    "    intpp2b = njit(lambda 𝜖, fk, 𝜃2, 𝜓2, 𝜉2, b: (w21(𝜖) + 𝜃2*(Y(𝜖, fk)-b) + 𝜉2)**(-𝜓2)*g(𝜖))\n",
    "    intqq2 = njit(lambda 𝜖, fk, 𝜃2, 𝜓2, b: (w21(𝜖) + 𝜃2*(Y(𝜖, fk)-b) + b)**(-𝜓2)*(Y(𝜖, fk) - b)*g(𝜖))\n",
    "\n",
    "\n",
    "    # Loop: Find fixed points V, q and p\n",
    "    V_crit = 1\n",
    "    while V_crit>1e-5:\n",
    "\n",
    "        # We begin by adding the guess for the value of the firm to endowment\n",
    "        V = (Vl+Vh)/2\n",
    "        ww10 = w10 + 𝜃10*V\n",
    "        ww20 = w20 + 𝜃20*V\n",
    "\n",
    "        # Production\n",
    "        fk = A*(k**𝛼)\n",
    "#         Y = lambda 𝜖: np.exp(𝜖)*fk\n",
    "\n",
    "        # Compute integration threshold\n",
    "        epstar = np.log(b/fk)\n",
    "\n",
    "\n",
    "        #**************************************************************\n",
    "        # Compute the prices and allocations consistent with consumers'\n",
    "        # Euler equations\n",
    "        #**************************************************************\n",
    "\n",
    "        # We impose the following:\n",
    "        # Agent 1 buys equity\n",
    "        # Agent 2 buys equity and all debt\n",
    "        # Agents trade such that prices converge\n",
    "\n",
    "        #========\n",
    "        # Agent 1\n",
    "        #========\n",
    "        # Holdings\n",
    "        𝜉1a = 0\n",
    "        𝜉1b = b/2\n",
    "        p = 0.3\n",
    "\n",
    "        while abs(𝜉1b - 𝜉1a) > 0.001:\n",
    "\n",
    "            𝜉1 = (𝜉1a + 𝜉1b) / 2\n",
    "            𝜃1a = 0.3\n",
    "            𝜃1b = 1\n",
    "\n",
    "            while abs(𝜃1b - 𝜃1a) > (0.001/b):\n",
    "\n",
    "                𝜃1 = (𝜃1a + 𝜃1b) / 2\n",
    "\n",
    "                # qq1 is the equity price consistent with agent-1 Euler Equation\n",
    "                ## Note: Price is in the date-0 budget constraint of the agent\n",
    "\n",
    "                ## First, compute the constant term that is not influenced by q\n",
    "                ## that is, 𝛽E[u'(c^{1}_{1})d^{e}(k,B)]\n",
    "#                 intqq1 = lambda 𝜖: (w11(𝜖) + 𝜃1*(Y(𝜖, fk) - b) + 𝜉1)**(-𝜓1)*(Y(𝜖, fk) - b)*g(𝜖)\n",
    "#                 const_qq1 = 𝛽 * quad(intqq1,epstar,bound)[0]\n",
    "                const_qq1 = 𝛽 * quad(intqq1,epstar,bound, args=(fk, 𝜃1, 𝜓1, 𝜉1, b))[0]\n",
    "\n",
    "                ## Second, iterate to get the equity price q\n",
    "                qq1l = 0\n",
    "                qq1h = ww10\n",
    "                diff = 1\n",
    "                while diff > 1e-7:\n",
    "                    qq1 = (qq1l+qq1h)/2\n",
    "                    rhs = const_qq1/((ww10-qq1*𝜃1-p*𝜉1)**(-𝜓1));\n",
    "                    if (rhs > qq1):\n",
    "                        qq1l = qq1\n",
    "                    else:\n",
    "                        qq1h = qq1\n",
    "                    diff = abs(qq1l-qq1h)\n",
    "\n",
    "                # pp1 is the bond price consistent with agent-2 Euler Equation\n",
    "                ## Note: Price is in the date-0 budget constraint of the agent\n",
    "\n",
    "                ## First, compute the constant term that is not influenced by p\n",
    "                ## that is, 𝛽E[u'(c^{1}_{1})d^{b}(k,B)]\n",
    "#                 intpp1a = lambda 𝜖: (Y(𝜖, fk)/b)*(w11(𝜖) + Y(𝜖, fk)/b*𝜉1)**(-𝜓1)*g(𝜖)\n",
    "#                 intpp1b = lambda 𝜖: (w11(𝜖) + 𝜃1*(Y(𝜖, fk)-b) + 𝜉1)**(-𝜓1)*g(𝜖)\n",
    "#                 const_pp1 = 𝛽 * (quad(intpp1a,-bound,epstar)[0] + quad(intpp1b,epstar,bound)[0])\n",
    "                const_pp1 = 𝛽 * (quad(intpp1a,-bound,epstar, args=(fk, 𝜓1, 𝜉1, b))[0] \\\n",
    "                                  + quad(intpp1b,epstar,bound, args=(fk, 𝜃1, 𝜓1, 𝜉1, b))[0])\n",
    "\n",
    "                ## iterate to get the bond price p\n",
    "                pp1l = 0\n",
    "                pp1h = ww10/b\n",
    "                diff = 1\n",
    "                while diff > 1e-7:\n",
    "                    pp1 = (pp1l+pp1h)/2\n",
    "                    rhs = const_pp1/((ww10-qq1*𝜃1-pp1*𝜉1)**(-𝜓1))\n",
    "                    if (rhs > pp1):\n",
    "                        pp1l = pp1\n",
    "                    else:\n",
    "                        pp1h = pp1\n",
    "                    diff = abs(pp1l-pp1h)\n",
    "\n",
    "                #========\n",
    "                # Agent 2\n",
    "                #========\n",
    "                𝜉2 = b - 𝜉1\n",
    "                𝜃2 = 1 - 𝜃1\n",
    "\n",
    "                # pp2 is the bond price consistent with agent-2 Euler Equation\n",
    "                ## Note: Price is in the date-0 budget constraint of the agent\n",
    "\n",
    "                ## First, compute the constant term that is not influenced by p\n",
    "                ## that is, 𝛽E[u'(c^{2}_{1})d^{b}(k,B)]\n",
    "#                 intpp2a = lambda 𝜖: (Y(𝜖, fk)/b)*(w21(𝜖) + Y(𝜖, fk)/b*𝜉2)**(-𝜓2)*g(𝜖)\n",
    "#                 intpp2b = lambda 𝜖: (w21(𝜖) + 𝜃2*(Y(𝜖, fk)-b) + 𝜉2)**(-𝜓2)*g(𝜖)\n",
    "#                 const_pp2 = 𝛽 * (quad(intpp2a,-bound,epstar)[0] + quad(intpp2b,epstar,bound)[0])\n",
    "                const_pp2 = 𝛽 * (quad(intpp2a,-bound,epstar, args=(fk, 𝜓2, 𝜉2, b))[0] \\\n",
    "                                  + quad(intpp2b,epstar,bound, args=(fk, 𝜃2, 𝜓2, 𝜉2, b))[0])\n",
    "\n",
    "                ## iterate to get the bond price p\n",
    "                pp2l = 0\n",
    "                pp2h = ww20/b\n",
    "                diff = 1\n",
    "                while diff > 1e-7:\n",
    "                    pp2 = (pp2l+pp2h)/2\n",
    "                    rhs = const_pp2/((ww20-qq1*𝜃2-pp2*𝜉2)**(-𝜓2))\n",
    "                    if (rhs > pp2):\n",
    "                        pp2l = pp2\n",
    "                    else:\n",
    "                        pp2h = pp2\n",
    "                    diff = abs(pp2l-pp2h)\n",
    "\n",
    "                # p be the maximum valuation for the bond among agents\n",
    "                ## This will be the equity price based on Makowski's criterion\n",
    "                p = max(pp1,pp2)\n",
    "\n",
    "\n",
    "                # qq2 is the equity price consistent with agent-2 Euler Equation\n",
    "#                 intqq2 = lambda 𝜖: (w21(𝜖) + 𝜃2*(Y(𝜖, fk)-b) + b)**(-𝜓2)*(Y(𝜖, fk) - b)*g(𝜖)\n",
    "#                 const_qq2 = 𝛽 * quad(intqq2,epstar,bound)[0]\n",
    "                const_qq2 = 𝛽 * quad(intqq2,epstar,bound, args=(fk, 𝜃2, 𝜓2, b))[0]\n",
    "                qq2l = 0\n",
    "                qq2h = ww20\n",
    "                diff = 1\n",
    "                while diff > 1e-7:\n",
    "                    qq2 = (qq2l+qq2h)/2\n",
    "                    rhs = const_qq2/((ww20-qq2*𝜃2-p*𝜉2)**(-𝜓2));\n",
    "                    if (rhs > qq2):\n",
    "                        qq2l = qq2\n",
    "                    else:\n",
    "                        qq2h = qq2\n",
    "                    diff = abs(qq2l-qq2h)\n",
    "\n",
    "                # q be the maximum valuation for the equity among agents\n",
    "                ## This will be the equity price based on Makowski's criterion\n",
    "                q = max(qq1,qq2)\n",
    "\n",
    "                #================\n",
    "                # Update holdings\n",
    "                #================\n",
    "                if qq1 > qq2:\n",
    "                    𝜃1a = 𝜃1\n",
    "                else:\n",
    "                    𝜃1b = 𝜃1\n",
    "\n",
    "                #print(p,q,𝜉1,𝜃1)\n",
    "\n",
    "            if pp1 > pp2:\n",
    "                𝜉1a = 𝜉1\n",
    "            else:\n",
    "                𝜉1b = 𝜉1\n",
    "\n",
    "\n",
    "        #================\n",
    "        # Get consumption\n",
    "        #================\n",
    "        c10 = ww10 - q*𝜃1 - p*𝜉1\n",
    "        c11 = lambda 𝜖: w11(𝜖) + 𝜃1*max(Y(𝜖, fk)-b,0) + 𝜉1*min(Y(𝜖, fk)/b,1)\n",
    "        c20 = ww20 - q*(1-𝜃1) - p*(b-𝜉1)\n",
    "        c21 = lambda 𝜖: w21(𝜖) + (1-𝜃1)*max(Y(𝜖, fk)-b,0) + (b-𝜉1)*min(Y(𝜖, fk)/b,1)\n",
    "\n",
    "        # Compute the value of the firm\n",
    "        value_x = -k + q + p*b\n",
    "        if (value_x > V):\n",
    "            Vl = V\n",
    "        else:\n",
    "            Vh = V\n",
    "        V_crit = abs(value_x-V)\n",
    "\n",
    "    return V,k,b,p,q,c10,c11,c20,c21,𝜉1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee0bd093",
   "metadata": {},
   "source": [
    "Here is our strategy for checking *stability* of an equilibrium.\n",
    "\n",
    "We use `off_eq_check` to obtain  consumption plans for both agents\n",
    "at the conjectured big $ K $ and big $ B $.\n",
    "\n",
    "Then we input  consumption plans into the function `eq_valuation`\n",
    "from the BCG model class and plot the agents’ valuations associated\n",
    "with different choices of $ k $ and $ b $.\n",
    "\n",
    "Our hunch is that $ (k^*,b^{**}) $ is **not** at the top of the\n",
    "firm valuation 3D surface so that the firm is **not** maximizing its\n",
    "value if it chooses $ k = K = k^* $ and $ b = B = b^{**} $.\n",
    "\n",
    "That indicates that $ (k^*,b^{**}) $ is not an equilibrium capital\n",
    "structure for the firm.\n",
    "\n",
    "We first check the case in which $ b^{**} = b^* - e $ where\n",
    "$ e = 0.1 $:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67752422",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "#====================== Experiment 1 ======================#\n",
    "Ve1,ke1,be1,pe1,qe1,c10e1,c11e1,c20e1,c21e1,𝜉1e1 = off_eq_check(mdl,\n",
    "                                                                kss,\n",
    "                                                                bss,\n",
    "                                                                e=-0.1)\n",
    "\n",
    "# Firm Valuation\n",
    "kgride1, bgride1, Vgride1, Qgride1, Pgride1 = mdl.eq_valuation(c10e1, c11e1, c20e1, c21e1,N=20)\n",
    "\n",
    "print('Maximum valuation of the firm value in the (k,b) grid: {:.4f}'.format(Vgride1.max()))\n",
    "print('Equilibrium firm value: {:.4f}'.format(Ve1))\n",
    "\n",
    "fig = go.Figure(data=[go.Scatter3d(x=[ke1],\n",
    "                                   y=[be1],\n",
    "                                   z=[Ve1],\n",
    "                                   mode='markers',\n",
    "                                   marker=dict(size=3, color='red')),\n",
    "                      go.Surface(x=kgride1,\n",
    "                                 y=bgride1,\n",
    "                                 z=Vgride1,\n",
    "                                 colorscale='Greens',opacity=0.6)])\n",
    "\n",
    "fig.update_layout(scene = dict(\n",
    "                    xaxis_title='x - Capital k',\n",
    "                    yaxis_title='y - Debt b',\n",
    "                    zaxis_title='z - Firm Value V',\n",
    "                    aspectratio = dict(x=1,y=1,z=1)),\n",
    "                  width=700,\n",
    "                  height=700,\n",
    "                  margin=dict(l=50, r=50, b=65, t=90))\n",
    "fig.update_layout(scene_camera=dict(eye=dict(x=1.5, y=-1.5, z=2)))\n",
    "fig.update_layout(title='Equilibrium firm valuation for the grid of (k,b)')\n",
    "\n",
    "\n",
    "# Export to PNG file\n",
    "Image(fig.to_image(format=\"png\"))\n",
    "# fig.show() will provide interactive plot when running\n",
    "# code locally"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70bd800f",
   "metadata": {},
   "source": [
    "In the above 3D surface of prospective firm valuations, the perturbed\n",
    "choice $ (k^*,b^{*}-e) $, represented by the red dot, is not at the\n",
    "top.\n",
    "\n",
    "The firm could issue more debts and attain a higher firm valuation from\n",
    "the market.\n",
    "\n",
    "Therefore, $ (k^*,b^{*}-e) $ would not be an equilibrium.\n",
    "\n",
    "Next, we check for $ b^{**} = b^* + e $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e852a38",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "#====================== Experiment 2 ======================#\n",
    "Ve2,ke2,be2,pe2,qe2,c10e2,c11e2,c20e2,c21e2,𝜉1e2 = off_eq_check(mdl,\n",
    "                                                                kss,\n",
    "                                                                bss,\n",
    "                                                                e=0.1)\n",
    "\n",
    "# Firm Valuation\n",
    "kgride2, bgride2, Vgride2, Qgride2, Pgride2 = mdl.eq_valuation(c10e2, c11e2, c20e2, c21e2,N=20)\n",
    "\n",
    "print('Maximum valuation of the firm value in the (k,b) grid: {:.4f}'.format(Vgride2.max()))\n",
    "print('Equilibrium firm value: {:.4f}'.format(Ve2))\n",
    "\n",
    "fig = go.Figure(data=[go.Scatter3d(x=[ke2],\n",
    "                                   y=[be2],\n",
    "                                   z=[Ve2],\n",
    "                                   mode='markers',\n",
    "                                   marker=dict(size=3, color='red')),\n",
    "                      go.Surface(x=kgride2,\n",
    "                                 y=bgride2,\n",
    "                                 z=Vgride2,\n",
    "                                 colorscale='Greens',opacity=0.6)])\n",
    "\n",
    "fig.update_layout(scene = dict(\n",
    "                    xaxis_title='x - Capital k',\n",
    "                    yaxis_title='y - Debt b',\n",
    "                    zaxis_title='z - Firm Value V',\n",
    "                    aspectratio = dict(x=1,y=1,z=1)),\n",
    "                  width=700,\n",
    "                  height=700,\n",
    "                  margin=dict(l=50, r=50, b=65, t=90))\n",
    "fig.update_layout(scene_camera=dict(eye=dict(x=1.5, y=-1.5, z=2)))\n",
    "fig.update_layout(title='Equilibrium firm valuation for the grid of (k,b)')\n",
    "\n",
    "# Export to PNG file\n",
    "Image(fig.to_image(format=\"png\"))\n",
    "# fig.show() will provide interactive plot when running\n",
    "# code locally"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2d4a24f",
   "metadata": {},
   "source": [
    "In contrast to $ (k^*,b^* - e) $, the 3D surface for\n",
    "$ (k^*,b^*+e) $ now indicates that a firm would want o *decrease*\n",
    "its debt issuance to attain a higher valuation.\n",
    "\n",
    "That incentive to deviate means that $ (k^*,b^*+e) $ is not an\n",
    "equilibrium capital structure for the firm.\n",
    "\n",
    "Interestingly, if consumers were to anticipate that firms would\n",
    "over-issue debt, i.e. $ B > b^* $, then both types of consumer would\n",
    "want to hold corporate debt.\n",
    "\n",
    "For example,  $ \\xi^1 > 0 $:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66a55e30",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "print('Bond holdings of agent 1: {:.3f}'.format(𝜉1e2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "011b1c25",
   "metadata": {},
   "source": [
    "Our two *stability experiments* suggest that the equilibrium capital\n",
    "structure $ (k^*,b^*) $ is locally unique even though **at the\n",
    "equilibrium** an individual firm would be willing to deviate from the\n",
    "representative firms’ equilibrium  debt choice.\n",
    "\n",
    "These experiments thus refine our discussion of the *qualified*\n",
    "Modigliani-Miller theorem that prevails in this example economy."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7680ccc8",
   "metadata": {},
   "source": [
    "#### Equilibrium equity and bond price functions\n",
    "\n",
    "It is also interesting to look at the equilibrium price functions\n",
    "$ q(k,b) $ and $ p(k,b) $ faced by firms in our rational\n",
    "expectations equilibrium."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "023270da",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Equity Valuation\n",
    "fig = go.Figure(data=[go.Scatter3d(x=[kss],\n",
    "                                   y=[bss],\n",
    "                                   z=[qss],\n",
    "                                   mode='markers',\n",
    "                                   marker=dict(size=3, color='red')),\n",
    "                      go.Surface(x=kgrid,\n",
    "                                 y=bgrid,\n",
    "                                 z=Qgrid,\n",
    "                                 colorscale='Blues',opacity=0.6)])\n",
    "\n",
    "fig.update_layout(scene = dict(\n",
    "                    xaxis_title='x - Capital k',\n",
    "                    yaxis_title='y - Debt b',\n",
    "                    zaxis_title='z - Equity price q',\n",
    "                    aspectratio = dict(x=1,y=1,z=1)),\n",
    "                  width=700,\n",
    "                  height=700,\n",
    "                  margin=dict(l=50, r=50, b=65, t=90))\n",
    "fig.update_layout(scene_camera=dict(eye=dict(x=1.5, y=-1.5, z=2)))\n",
    "fig.update_layout(title='Equilibrium equity valuation for the grid of (k,b)')\n",
    "\n",
    "\n",
    "# Export to PNG file\n",
    "Image(fig.to_image(format=\"png\"))\n",
    "# fig.show() will provide interactive plot when running\n",
    "# code locally"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ce67911",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Bond Valuation\n",
    "fig = go.Figure(data=[go.Scatter3d(x=[kss],\n",
    "                                   y=[bss],\n",
    "                                   z=[pss],\n",
    "                                   mode='markers',\n",
    "                                   marker=dict(size=3, color='red')),\n",
    "                      go.Surface(x=kgrid,\n",
    "                                 y=bgrid,\n",
    "                                 z=Pgrid,\n",
    "                                 colorscale='Oranges',opacity=0.6)])\n",
    "\n",
    "fig.update_layout(scene = dict(\n",
    "                    xaxis_title='x - Capital k',\n",
    "                    yaxis_title='y - Debt b',\n",
    "                    zaxis_title='z - Bond price q',\n",
    "                    aspectratio = dict(x=1,y=1,z=1)),\n",
    "                  width=700,\n",
    "                  height=700,\n",
    "                  margin=dict(l=50, r=50, b=65, t=90))\n",
    "fig.update_layout(scene_camera=dict(eye=dict(x=1.5, y=-1.5, z=2)))\n",
    "fig.update_layout(title='Equilibrium bond valuation for the grid of (k,b)')\n",
    "\n",
    "\n",
    "# Export to PNG file\n",
    "Image(fig.to_image(format=\"png\"))\n",
    "# fig.show() will provide interactive plot when running\n",
    "# code locally"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aaeeda74",
   "metadata": {},
   "source": [
    "### Comments on equilibrium pricing functions\n",
    "\n",
    "The equilibrium pricing functions displayed above merit study and\n",
    "reflection.\n",
    "\n",
    "They reveal the countervailing effects on a firm’s valuations of bonds\n",
    "and equities that lie beneath the Modigliani-Miller ridge apparent in\n",
    "our earlier graph of an individual firm $ \\zeta $’s value as a\n",
    "function of $ k(\\zeta), b(\\zeta) $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8924083c",
   "metadata": {},
   "source": [
    "### Another example economy\n",
    "\n",
    "We illustrate how the fraction of initial endowments held by agent 2,\n",
    "$ w^2_0/(w^1_0+w^2_0) $ affects an equilibrium capital structure\n",
    "$ (k,b) = (K, B) $ well as associated equilibrium allocations.\n",
    "\n",
    "We are interested in how  agents 1 and 2\n",
    "value equity and  bond.\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "Q^i = \\beta \\int \\frac{u^\\prime(C^{i,*}_1(\\epsilon))}{u^\\prime(C^{i,*}_0)} d^e(k^*,b^*;\\epsilon) g(\\epsilon) \\ d\\epsilon \\\\\n",
    "P^i = \\beta \\int \\frac{u^\\prime(C^{i,*}_1(\\epsilon))}{u^\\prime(C^{i,*}_0)} d^b(k^*,b^*;\\epsilon) g(\\epsilon) \\ d\\epsilon \\\\\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "The function `valuations_by_agent` is used in calculating these\n",
    "valuations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea3a71f7",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Lists for storage\n",
    "wlist = []\n",
    "klist = []\n",
    "blist = []\n",
    "qlist = []\n",
    "plist = []\n",
    "Vlist = []\n",
    "tlist = []\n",
    "q1list = []\n",
    "q2list = []\n",
    "p1list = []\n",
    "p2list = []\n",
    "\n",
    "# For loop: optimization for each endowment combination\n",
    "for i in range(10):\n",
    "    print(i)\n",
    "\n",
    "    # Save fraction\n",
    "    w10 = 0.9 - 0.05*i\n",
    "    w20 = 1.1 + 0.05*i\n",
    "    wlist.append(w20/(w10+w20))\n",
    "\n",
    "    # Create the instance\n",
    "    mdl = BCG_incomplete_markets(w10 = w10, w20 = w20, ktop = 0.5, btop = 2.5)\n",
    "\n",
    "    # Solve for equilibrium\n",
    "    kss,bss,Vss,qss,pss,c10ss,c11ss,c20ss,c21ss,𝜃1ss = mdl.solve_eq(print_crit=False)\n",
    "\n",
    "    # Store the equilibrium results\n",
    "    klist.append(kss)\n",
    "    blist.append(bss)\n",
    "    qlist.append(qss)\n",
    "    plist.append(pss)\n",
    "    Vlist.append(Vss)\n",
    "    tlist.append(𝜃1ss)\n",
    "\n",
    "    # Evaluations of equity and bond by each agent\n",
    "    Q1,Q2,P1,P2 = mdl.valuations_by_agent(c10ss, c11ss, c20ss, c21ss, kss, bss)\n",
    "\n",
    "    # Save the valuations\n",
    "    q1list.append(Q1)\n",
    "    q2list.append(Q2)\n",
    "    p1list.append(P1)\n",
    "    p2list.append(P2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c5a5345",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Plot\n",
    "fig, ax = plt.subplots(3,2,figsize=(12,12))\n",
    "ax[0,0].plot(wlist,klist)\n",
    "ax[0,0].set_title('capital')\n",
    "ax[0,1].plot(wlist,blist)\n",
    "ax[0,1].set_title('debt')\n",
    "ax[1,0].plot(wlist,qlist)\n",
    "ax[1,0].set_title('equity price')\n",
    "ax[1,1].plot(wlist,plist)\n",
    "ax[1,1].set_title('bond price')\n",
    "ax[2,0].plot(wlist,Vlist)\n",
    "ax[2,0].set_title('firm value')\n",
    "ax[2,0].set_xlabel('fraction of initial endowment held by agent 2',fontsize=13)\n",
    "\n",
    "# Create a list of Default thresholds\n",
    "A = mdl.A\n",
    "𝛼 = mdl.𝛼\n",
    "epslist = []\n",
    "for i in range(len(wlist)):\n",
    "    bb = blist[i]\n",
    "    kk = klist[i]\n",
    "    eps = np.log(bb/(A*kk**𝛼))\n",
    "    epslist.append(eps)\n",
    "\n",
    "# Plot (cont.)\n",
    "ax[2,1].plot(wlist,epslist)\n",
    "ax[2,1].set_title(r'default threshold $\\epsilon^*$')\n",
    "ax[2,1].set_xlabel('fraction of initial endowment held by agent 2',fontsize=13)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0af0b2d8",
   "metadata": {},
   "source": [
    "## A picture worth a thousand words\n",
    "\n",
    "Please stare at the above panels.\n",
    "\n",
    "They describe how equilibrium prices and quantities respond to\n",
    "alterations in the structure of society’s *hedging desires* across\n",
    "economies with different allocations of the initial endowment to our two\n",
    "types of agents.\n",
    "\n",
    "Now let’s see how the two types of agents value bonds and equities,\n",
    "keeping in mind that the type that values the asset highest determines\n",
    "the equilibrium price (and thus the pertinent set of Big $ C $’s)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26c3986f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Comparing the prices\n",
    "fig, ax = plt.subplots(1,3,figsize=(16,6))\n",
    "\n",
    "ax[0].plot(wlist,q1list,label='agent 1',color='green')\n",
    "ax[0].plot(wlist,q2list,label='agent 2',color='blue')\n",
    "ax[0].plot(wlist,qlist,label='equity price',color='red',linestyle='--')\n",
    "ax[0].legend()\n",
    "ax[0].set_title('equity valuations')\n",
    "ax[0].set_xlabel('fraction of initial endowment held by agent 2',fontsize=11)\n",
    "\n",
    "ax[1].plot(wlist,p1list,label='agent 1',color='green')\n",
    "ax[1].plot(wlist,p2list,label='agent 2',color='blue')\n",
    "ax[1].plot(wlist,plist,label='bond price',color='red',linestyle='--')\n",
    "ax[1].legend()\n",
    "ax[1].set_title('bond valuations')\n",
    "ax[1].set_xlabel('fraction of initial endowment held by agent 2',fontsize=11)\n",
    "\n",
    "ax[2].plot(wlist,tlist,color='blue')\n",
    "ax[2].set_title('equity holdings by agent 1')\n",
    "ax[2].set_xlabel('fraction of initial endowment held by agent 2',fontsize=11)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9af167c",
   "metadata": {},
   "source": [
    "It is rewarding to stare at the above plots too.\n",
    "\n",
    "In equilibrium, equity valuations are the same across the two types of\n",
    "agents but bond valuations are not.\n",
    "\n",
    "Agents of type 2 value bonds more highly (they want more hedging).\n",
    "\n",
    "Taken together with our earlier plot of equity holdings, these graphs confirm our earlier conjecture that while both type\n",
    "of agents hold equities, only agents of type 2 holds bonds."
   ]
  }
 ],
 "metadata": {
  "date": 1734843459.8564742,
  "filename": "BCG_incomplete_mkts.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Equilibrium Capital Structures with Incomplete Markets"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}