arellano.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<a id='arellano'></a>\n",
    "<div id=\"qe-notebook-header\" align=\"right\" style=\"text-align:right;\">\n",
    "        <a href=\"https://quantecon.org/\" title=\"quantecon.org\">\n",
    "                <img style=\"width:250px;display:inline;\" width=\"250px\" src=\"https://assets.quantecon.org/img/qe-menubar-logo.svg\" alt=\"QuantEcon\">\n",
    "        </a>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Default Risk and Income Fluctuations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Default Risk and Income Fluctuations](#Default-Risk-and-Income-Fluctuations)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [Structure](#Structure)  \n",
    "  - [Equilibrium](#Equilibrium)  \n",
    "  - [Computation](#Computation)  \n",
    "  - [Results](#Results)  \n",
    "  - [Exercises](#Exercises)  \n",
    "  - [Solutions](#Solutions)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition to what&rsquo;s in Anaconda, this lecture will need the following libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": true
   },
   "outputs": [],
   "source": [
    "!pip install --upgrade quantecon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "This lecture computes versions of  Arellano&rsquo;s  [[Are08]](https://lectures.quantecon.org/zreferences.html#arellano2008default) model of sovereign default.\n",
    "\n",
    "The model describes interactions among default risk, output,  and  an equilibrium interest rate  that includes a premium for endogenous default risk.\n",
    "\n",
    "The decision maker is a government of a small open economy that borrows from risk-neutral foreign creditors.\n",
    "\n",
    "The foreign lenders must be compensated for default risk.\n",
    "\n",
    "The government borrows and lends abroad in order to  smooth the consumption of its citizens.\n",
    "\n",
    "The government repays its debt only if it wants to, but declining to pay has adverse consequences.\n",
    "\n",
    "The interest rate on government debt adjusts in response to the state-dependent default probability chosen by government.\n",
    "\n",
    "The model yields outcomes that help interpret sovereign default experiences, including\n",
    "\n",
    "- countercyclical interest rates on sovereign debt  \n",
    "- countercyclical trade balances  \n",
    "- high volatility of consumption relative to output  \n",
    "\n",
    "\n",
    "Notably, long recessions caused by bad draws in the income process increase the government&rsquo;s incentive to default.\n",
    "\n",
    "This can lead to\n",
    "\n",
    "- spikes in interest rates  \n",
    "- temporary losses of access to international credit markets  \n",
    "- large drops in output, consumption, and welfare  \n",
    "- large capital outflows during recessions  \n",
    "\n",
    "\n",
    "Such dynamics are consistent with experiences of many countries.\n",
    "\n",
    "Let&rsquo;s start with some imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import quantecon as qe\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from numba import jit\n",
    "import random"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Structure\n",
    "\n",
    "In this section we describe the main features of the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Output, Consumption and Debt\n",
    "\n",
    "A small open economy is endowed with an exogenous stochastically fluctuating potential output stream $ \\{y_t\\} $.\n",
    "\n",
    "Potential output is realized only in periods in which the government honors its sovereign debt.\n",
    "\n",
    "The output good can be traded or consumed.\n",
    "\n",
    "The sequence $ \\{y_t\\} $ is described by a Markov process with stochastic density kernel $ p(y, y') $.\n",
    "\n",
    "Households within the country are identical and rank stochastic consumption streams according to\n",
    "\n",
    "\n",
    "<a id='equation-utility'></a>\n",
    "$$\n",
    "\\mathbb E \\sum_{t=0}^{\\infty} \\beta^t u(c_t) \\tag{1}\n",
    "$$\n",
    "\n",
    "Here\n",
    "\n",
    "- $ 0 < \\beta < 1 $ is a time discount factor  \n",
    "- $ u $ is an increasing and strictly concave utility function  \n",
    "\n",
    "\n",
    "Consumption sequences enjoyed by households are affected by the government&rsquo;s decision to borrow or lend internationally.\n",
    "\n",
    "The government is benevolent in the sense that its aim is to maximize [(1)](#equation-utility).\n",
    "\n",
    "The government is the only domestic actor with access to foreign credit.\n",
    "\n",
    "Because household are averse to consumption fluctuations, the government will try to smooth consumption by borrowing from (and lending to) foreign creditors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Asset Markets\n",
    "\n",
    "The only credit instrument available to the government is a one-period bond traded in international credit markets.\n",
    "\n",
    "The bond market has the following features\n",
    "\n",
    "- The bond matures in one period and is not state contingent.  \n",
    "- A purchase of a bond with face value $ B' $ is a claim to $ B' $ units of the\n",
    "  consumption good next period.  \n",
    "- To purchase $ B' $  next period costs $ q B' $ now, or, what is equivalent.  \n",
    "- For selling $ -B' $ units of next period goods the seller earns $ - q B' $ of today&rsquo;s goods.  \n",
    "  \n",
    "  - If $ B' < 0 $, then $ -q B' $ units of the good are received in the current period, for a promise to repay $ -B' $ units next period.  \n",
    "  - There is an equilibrium  price function $ q(B', y) $ that makes $ q $ depend on both $ B' $ and $ y $.  \n",
    "  \n",
    "\n",
    "\n",
    "Earnings on the government portfolio are distributed (or, if negative, taxed) lump sum to households.\n",
    "\n",
    "When the government is not excluded from financial markets, the one-period national budget constraint is\n",
    "\n",
    "\n",
    "<a id='equation-resource'></a>\n",
    "$$\n",
    "c = y + B - q(B', y) B' \\tag{2}\n",
    "$$\n",
    "\n",
    "Here and below, a prime denotes a next period value or a claim maturing next period.\n",
    "\n",
    "To rule out Ponzi schemes, we also require that $ B \\geq -Z $ in every period.\n",
    "\n",
    "- $ Z $ is chosen to be sufficiently large that the constraint never binds in equilibrium.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Financial Markets\n",
    "\n",
    "Foreign creditors\n",
    "\n",
    "- are risk neutral  \n",
    "- know the domestic output stochastic process $ \\{y_t\\} $ and observe $ y_t, y_{t-1}, \\ldots, $ at time $ t $  \n",
    "- can borrow or lend without limit in an international credit market at a constant international interest rate $ r $  \n",
    "- receive full payment if the government chooses to pay  \n",
    "- receive zero if the government defaults on its one-period debt due  \n",
    "\n",
    "\n",
    "When a government is expected to default next period with  probability $ \\delta $,  the expected value of a promise to pay one unit of consumption next period is $ 1 - \\delta $.\n",
    "\n",
    "Therefore, the discounted expected value of a promise to pay $ B $ next period is\n",
    "\n",
    "\n",
    "<a id='equation-epc'></a>\n",
    "$$\n",
    "q = \\frac{1 - \\delta}{1 + r} \\tag{3}\n",
    "$$\n",
    "\n",
    "Next we turn to how the government in effect chooses the default probability $ \\delta $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Government&rsquo;s Decisions\n",
    "\n",
    "At each point in time $ t $, the government chooses between\n",
    "\n",
    "1. defaulting  \n",
    "1. meeting its current obligations and purchasing or selling an optimal quantity of  one-period sovereign debt  \n",
    "\n",
    "\n",
    "Defaulting means declining to repay all of  its current obligations.\n",
    "\n",
    "If the government defaults in the current period, then consumption equals current output.\n",
    "\n",
    "But a sovereign default has two consequences:\n",
    "\n",
    "1. Output immediately falls from $ y $ to $ h(y) $, where $ 0 \\leq h(y) \\leq y $.  \n",
    "  \n",
    "  - It returns to $ y $ only after the country regains access to  international credit markets.  \n",
    "  \n",
    "1. The country loses access to foreign credit markets.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reentering International Credit Market\n",
    "\n",
    "While in a state of default, the economy regains access to\n",
    "foreign credit in each subsequent period with probability\n",
    "$ \\theta $."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equilibrium\n",
    "\n",
    "Informally, an equilibrium is a sequence of interest rates on its sovereign debt, a stochastic sequence of government default decisions  and an implied flow of household consumption such that\n",
    "\n",
    "1. Consumption and assets satisfy the national budget constraint.  \n",
    "1. The government maximizes household utility taking into account  \n",
    "  - the resource constraint  \n",
    "  - the effect of its choices on the price of bonds  \n",
    "  - consequences of defaulting now for future net output and future borrowing and lending opportunities  \n",
    "1. The interest rate on the government&rsquo;s debt includes a risk-premium sufficient to make foreign creditors expect on average to earn the constant risk-free international interest rate.  \n",
    "\n",
    "\n",
    "To express these ideas more precisely, consider first the choices of the\n",
    "government, which\n",
    "\n",
    "1. enters a period with initial assets $ B $, or  what is the same thing, initial debt to be repaid now of $ -B $  \n",
    "1. observes current output $ y $, and  \n",
    "1. chooses either  \n",
    "  \n",
    "  1. to default, or  \n",
    "  1. to pay  $ -B $ and set next period&rsquo;s debt due to $ -B' $  \n",
    "  \n",
    "\n",
    "\n",
    "In a  recursive formulation,\n",
    "\n",
    "- state variables for the government comprise the pair $ (B, y) $  \n",
    "- $ v(B, y) $ is the optimum value of the government&rsquo;s problem when at the beginning of a period  it faces the choice of whether to honor or default  \n",
    "- $ v_c(B, y) $ is the value of choosing to pay obligations falling due  \n",
    "- $ v_d(y) $ is the value of choosing to default  \n",
    "\n",
    "\n",
    "$ v_d(y) $ does not depend on $ B $ because, when access to credit is eventually regained, net foreign assets equal $ 0 $.\n",
    "\n",
    "Expressed recursively, the value of defaulting is\n",
    "\n",
    "$$\n",
    "v_d(y) = u(h(y)) +\n",
    "            \\beta \\int \\left\\{\n",
    "            \\theta v(0, y') + (1 - \\theta) v_d(y')\n",
    "            \\right\\}\n",
    "            p(y, y') dy'\n",
    "$$\n",
    "\n",
    "The value of paying is\n",
    "\n",
    "$$\n",
    "v_c(B, y) = \\max_{B' \\geq -Z}\n",
    "       \\left\\{\n",
    "            u(y - q(B', y) B' + B) +\n",
    "            \\beta \\int v(B', y') p(y, y') dy'\n",
    "      \\right\\}\n",
    "$$\n",
    "\n",
    "The three value functions are linked by\n",
    "\n",
    "$$\n",
    "v(B, y) = \\max\\{ v_c(B, y), v_d(y) \\}\n",
    "$$\n",
    "\n",
    "The government chooses to default when\n",
    "\n",
    "$$\n",
    "v_c(B, y) < v_d(y)\n",
    "$$\n",
    "\n",
    "and hence given $ B' $ the probability of default next period  is\n",
    "\n",
    "\n",
    "<a id='equation-delta'></a>\n",
    "$$\n",
    "\\delta(B', y) := \\int \\mathbb 1\\{v_c(B', y') < v_d(y') \\} p(y, y') dy' \\tag{4}\n",
    "$$\n",
    "\n",
    "Given zero profits for foreign creditors in equilibrium, we can combine [(3)](#equation-epc) and [(4)](#equation-delta) to pin down\n",
    "the bond price function:\n",
    "\n",
    "\n",
    "<a id='equation-bondprice'></a>\n",
    "$$\n",
    "q(B', y) = \\frac{1 - \\delta(B', y)}{1 + r} \\tag{5}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition of Equilibrium\n",
    "\n",
    "An *equilibrium* is\n",
    "\n",
    "- a  pricing function $ q(B',y) $,  \n",
    "- a triple of value functions $ (v_c(B, y), v_d(y), v(B,y)) $,  \n",
    "- a decision rule telling the government when to default and when to pay as a function of the state  $ (B, y) $, and  \n",
    "- an asset accumulation rule that, conditional on choosing not to  default,  maps $ (B,y) $ into $ B' $  \n",
    "\n",
    "\n",
    "such that\n",
    "\n",
    "- The three Bellman equations for $ (v_c(B, y), v_d(y), v(B,y)) $ are satisfied  \n",
    "- Given the price function $ q(B',y) $, the default decision rule and the asset accumulation decision rule attain the optimal value function  $ v(B,y) $, and  \n",
    "- The price function $ q(B',y) $ satisfies equation [(5)](#equation-bondprice)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computation\n",
    "\n",
    "Let&rsquo;s now compute an equilibrium of Arellano&rsquo;s model.\n",
    "\n",
    "The equilibrium objects are the value function $ v(B, y) $, the associated default decision rule, and the pricing function $ q(B', y) $.\n",
    "\n",
    "We&rsquo;ll use our code to replicate Arellano&rsquo;s results.\n",
    "\n",
    "After that we&rsquo;ll perform some additional simulations.\n",
    "\n",
    "The majority of the code below was written by [Chase Coleman](https://github.com/cc7768).\n",
    "\n",
    "It uses a slightly modified version of the algorithm recommended by Arellano.\n",
    "\n",
    "- The appendix to [[Are08]](https://lectures.quantecon.org/zreferences.html#arellano2008default) recommends\n",
    "  value function iteration until convergence, updating the price, and then repeating.  \n",
    "- Instead, we update the bond price at every value function iteration step.  \n",
    "\n",
    "\n",
    "The second approach is faster and the two different procedures deliver very similar results.\n",
    "\n",
    "Here is a more detailed description of our algorithm:\n",
    "\n",
    "1. Guess a value function $ v(B, y) $ and price function $ q(B', y) $.  \n",
    "1. At each pair $ (B, y) $,  \n",
    "  - update the value of defaulting $ v_d(y) $.  \n",
    "  - update the value of continuing $ v_c(B, y) $.  \n",
    "1. Update the value function $ v(B, y) $, the default rule,  the implied ex ante default probability,  and the price function.  \n",
    "1. Check for convergence. If converged, stop &ndash; if not, go to step 2.  \n",
    "\n",
    "\n",
    "We use simple discretization on a grid of asset holdings and income levels.\n",
    "\n",
    "The output process is discretized using [Tauchen&rsquo;s quadrature method](https://github.com/QuantEcon/QuantEcon.py/blob/master/quantecon/markov/approximation.py).\n",
    "\n",
    "[Numba](https://lectures.quantecon.org/numba.html#numba-link) has been used in two places to speed up the code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "\n",
    "Authors: Chase Coleman, John Stachurski\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "class Arellano_Economy:\n",
    "    \"\"\"\n",
    "    Arellano 2008 deals with a small open economy whose government\n",
    "    invests in foreign assets in order to smooth the consumption of\n",
    "    domestic households. Domestic households receive a stochastic\n",
    "    path of income.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    &beta; : float\n",
    "        Time discounting parameter\n",
    "    &gamma; : float\n",
    "        Risk-aversion parameter\n",
    "    r : float\n",
    "        int lending rate\n",
    "    &rho; : float\n",
    "        Persistence in the income process\n",
    "    &eta; : float\n",
    "        Standard deviation of the income process\n",
    "    &theta; : float\n",
    "        Probability of re-entering financial markets in each period\n",
    "    ny : int\n",
    "        Number of points in y grid\n",
    "    nB : int\n",
    "        Number of points in B grid\n",
    "    tol : float\n",
    "        Error tolerance in iteration\n",
    "    maxit : int\n",
    "        Maximum number of iterations\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self,\n",
    "                &beta;=.953,         # time discount rate\n",
    "                &gamma;=2.,           # risk aversion\n",
    "                r=0.017,        # international interest rate\n",
    "                &rho;=.945,         # persistence in output\n",
    "                &eta;=0.025,        # st dev of output shock\n",
    "                &theta;=0.282,        # prob of regaining access\n",
    "                ny=21,          # number of points in y grid\n",
    "                nB=251,         # number of points in B grid\n",
    "                tol=1e-8,       # error tolerance in iteration\n",
    "                maxit=10000):\n",
    "\n",
    "        # Save parameters\n",
    "        self.&beta;, self.&gamma;, self.r = &beta;, &gamma;, r\n",
    "        self.&rho;, self.&eta;, self.&theta; = &rho;, &eta;, &theta;\n",
    "        self.ny, self.nB = ny, nB\n",
    "\n",
    "        # Create grids and discretize Markov process\n",
    "        self.Bgrid = np.linspace(-.45, .45, nB)\n",
    "        self.mc = qe.markov.tauchen(&rho;, &eta;, 0, 3, ny)\n",
    "        self.ygrid = np.exp(self.mc.state_values)\n",
    "        self.Py = self.mc.P\n",
    "\n",
    "        # Output when in default\n",
    "        ymean = np.mean(self.ygrid)\n",
    "        self.def_y = np.minimum(0.969 * ymean, self.ygrid)\n",
    "\n",
    "        # Allocate memory\n",
    "        self.Vd = np.zeros(ny)\n",
    "        self.Vc = np.zeros((ny, nB))\n",
    "        self.V = np.zeros((ny, nB))\n",
    "        self.Q = np.ones((ny, nB)) * .95  # Initial guess for prices\n",
    "        self.default_prob = np.empty((ny, nB))\n",
    "\n",
    "        # Compute the value functions, prices, and default prob\n",
    "        self.solve(tol=tol, maxit=maxit)\n",
    "        # Compute the optimal savings policy conditional on no default\n",
    "        self.compute_savings_policy()\n",
    "\n",
    "    def solve(self, tol=1e-8, maxit=10000):\n",
    "        # Iteration Stuff\n",
    "        it = 0\n",
    "        dist = 10.\n",
    "\n",
    "        # Alloc memory to store next iterate of value function\n",
    "        V_upd = np.zeros((self.ny, self.nB))\n",
    "\n",
    "        # Main loop\n",
    "        while dist > tol and maxit > it:\n",
    "\n",
    "            # Compute expectations for this iteration\n",
    "            Vs = self.V, self.Vd, self.Vc\n",
    "            EV, EVd, EVc = (self.Py @ v for v in Vs)\n",
    "\n",
    "            # Run inner loop to update value functions Vc and Vd.\n",
    "            # Note that Vc and Vd are updated in place.  Other objects\n",
    "            # are not modified.\n",
    "            _inner_loop(self.ygrid, self.def_y,\n",
    "                        self.Bgrid, self.Vd, self.Vc,\n",
    "                        EVc, EVd, EV, self.Q,\n",
    "                        self.&beta;, self.&theta;, self.&gamma;)\n",
    "\n",
    "            # Update prices\n",
    "            Vd_compat = np.repeat(self.Vd, self.nB).reshape(self.ny, self.nB)\n",
    "            default_states = Vd_compat > self.Vc\n",
    "            self.default_prob[:, :] = self.Py @ default_states\n",
    "            self.Q[:, :] = (1 - self.default_prob)/(1 + self.r)\n",
    "\n",
    "            # Update main value function and distance\n",
    "            V_upd[:, :] = np.maximum(self.Vc, Vd_compat)\n",
    "            dist = np.max(np.abs(V_upd - self.V))\n",
    "            self.V[:, :] = V_upd[:, :]\n",
    "\n",
    "            it += 1\n",
    "            if it % 25 == 0:\n",
    "                print(f\"Running iteration {it} with dist of {dist}\")\n",
    "\n",
    "        return None\n",
    "\n",
    "    def compute_savings_policy(self):\n",
    "        \"\"\"\n",
    "        Compute optimal savings B' conditional on not defaulting.\n",
    "        The policy is recorded as an index value in Bgrid.\n",
    "        \"\"\"\n",
    "\n",
    "        # Allocate memory\n",
    "        self.next_B_index = np.empty((self.ny, self.nB))\n",
    "        EV = self.Py @ self.V\n",
    "\n",
    "        _compute_savings_policy(self.ygrid, self.Bgrid, self.Q, EV,\n",
    "                                self.&gamma;, self.&beta;, self.next_B_index)\n",
    "\n",
    "    def simulate(self, T, y_init=None, B_init=None):\n",
    "        \"\"\"\n",
    "        Simulate time series for output, consumption, B'.\n",
    "        \"\"\"\n",
    "        # Find index i such that Bgrid[i] is near 0\n",
    "        zero_B_index = np.searchsorted(self.Bgrid, 0)\n",
    "\n",
    "        if y_init is None:\n",
    "            # Set to index near the mean of the ygrid\n",
    "            y_init = np.searchsorted(self.ygrid, self.ygrid.mean())\n",
    "        if B_init is None:\n",
    "            B_init = zero_B_index\n",
    "        # Start off not in default\n",
    "        in_default = False\n",
    "\n",
    "        y_sim_indices = self.mc.simulate_indices(T, init=y_init)\n",
    "        B_sim_indices = np.empty(T, dtype=np.int64)\n",
    "        B_sim_indices[0] = B_init\n",
    "        q_sim = np.empty(T)\n",
    "        in_default_series = np.zeros(T, dtype=np.int64)\n",
    "\n",
    "        for t in range(T-1):\n",
    "            yi, Bi = y_sim_indices[t], B_sim_indices[t]\n",
    "            if not in_default:\n",
    "                if self.Vc[yi, Bi] < self.Vd[yi]:\n",
    "                    in_default = True\n",
    "                    Bi_next = zero_B_index\n",
    "                else:\n",
    "                    new_index = self.next_B_index[yi, Bi]\n",
    "                    Bi_next = new_index\n",
    "            else:\n",
    "                in_default_series[t] = 1\n",
    "                Bi_next = zero_B_index\n",
    "                if random.uniform(0, 1) < self.&theta;:\n",
    "                    in_default = False\n",
    "            B_sim_indices[t+1] = Bi_next\n",
    "            q_sim[t] = self.Q[yi, int(Bi_next)]\n",
    "\n",
    "        q_sim[-1] = q_sim[-2]  # Extrapolate for the last price\n",
    "        return_vecs = (self.ygrid[y_sim_indices],\n",
    "                    self.Bgrid[B_sim_indices],\n",
    "                    q_sim,\n",
    "                    in_default_series)\n",
    "\n",
    "        return return_vecs\n",
    "\n",
    "\n",
    "@jit(nopython=True)\n",
    "def u(c, &gamma;):\n",
    "    return c**(1-&gamma;)/(1-&gamma;)\n",
    "\n",
    "\n",
    "@jit(nopython=True)\n",
    "def _inner_loop(ygrid, def_y, Bgrid, Vd, Vc, EVc,\n",
    "                EVd, EV, qq, &beta;, &theta;, &gamma;):\n",
    "    \"\"\"\n",
    "    This is a numba version of the inner loop of the solve in the\n",
    "    Arellano class. It updates Vd and Vc in place.\n",
    "    \"\"\"\n",
    "    ny, nB = len(ygrid), len(Bgrid)\n",
    "    zero_ind = nB // 2  # Integer division\n",
    "    for iy in range(ny):\n",
    "        y = ygrid[iy]   # Pull out current y\n",
    "\n",
    "        # Compute Vd\n",
    "        Vd[iy] = u(def_y[iy], &gamma;) + \\\n",
    "                &beta; * (&theta; * EVc[iy, zero_ind] + (1 - &theta;) * EVd[iy])\n",
    "\n",
    "        # Compute Vc\n",
    "        for ib in range(nB):\n",
    "            B = Bgrid[ib]  # Pull out current B\n",
    "\n",
    "            current_max = -1e14\n",
    "            for ib_next in range(nB):\n",
    "                c = max(y - qq[iy, ib_next] * Bgrid[ib_next] + B, 1e-14)\n",
    "                m = u(c, &gamma;) + &beta; * EV[iy, ib_next]\n",
    "                if m > current_max:\n",
    "                    current_max = m\n",
    "            Vc[iy, ib] = current_max\n",
    "\n",
    "    return None\n",
    "\n",
    "\n",
    "@jit(nopython=True)\n",
    "def _compute_savings_policy(ygrid, Bgrid, Q, EV, &gamma;, &beta;, next_B_index):\n",
    "    # Compute best index in Bgrid given iy, ib\n",
    "    ny, nB = len(ygrid), len(Bgrid)\n",
    "    for iy in range(ny):\n",
    "        y = ygrid[iy]\n",
    "        for ib in range(nB):\n",
    "            B = Bgrid[ib]\n",
    "            current_max = -1e10\n",
    "            for ib_next in range(nB):\n",
    "                c = max(y - Q[iy, ib_next] * Bgrid[ib_next] + B, 1e-14)\n",
    "                m = u(c, &gamma;) + &beta; * EV[iy, ib_next]\n",
    "                if m > current_max:\n",
    "                    current_max = m\n",
    "                    current_max_index = ib_next\n",
    "            next_B_index[iy, ib] = current_max_index\n",
    "    return None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results\n",
    "\n",
    "Let&rsquo;s start by trying to replicate the results obtained in [[Are08]](https://lectures.quantecon.org/zreferences.html#arellano2008default).\n",
    "\n",
    "In what follows, all results are computed using Arellano&rsquo;s parameter values.\n",
    "\n",
    "The values can be seen in the `__init__` method of the `Arellano_Economy` shown above.\n",
    "\n",
    "- For example, `r=0.017` matches the average quarterly rate on a 5 year US treasury\n",
    "  over the period 1983&ndash;2001.  \n",
    "\n",
    "\n",
    "Details on how to compute the figures are reported as solutions to the exercises.\n",
    "\n",
    "The first figure shows the bond price schedule and replicates Figure 3 of Arellano, where $ y_L $ and $ Y_H $ are particular below average and above average values of output $ y $.\n",
    "\n",
    "<img src=\"https://s3-ap-southeast-2.amazonaws.com/lectures.quantecon.org/py/_static/lecture_specific/arellano/arellano_bond_prices.png\" style=\"\">\n",
    "\n",
    "  \n",
    "- $ y_L $ is 5% below the mean of the $ y $ grid values  \n",
    "- $ y_H $ is 5% above  the mean of the $ y $ grid values  \n",
    "\n",
    "\n",
    "The grid used to compute this figure was relatively coarse (`ny, nB = 21, 251`) in order to match Arrelano&rsquo;s findings.\n",
    "\n",
    "Here&rsquo;s the same relationships computed on a finer grid (`ny, nB = 51, 551`)\n",
    "\n",
    "<img src=\"https://s3-ap-southeast-2.amazonaws.com/lectures.quantecon.org/py/_static/lecture_specific/arellano/arellano_bond_prices_2.png\" style=\"\">\n",
    "\n",
    "  \n",
    "In either case, the figure shows that\n",
    "\n",
    "- Higher levels of debt (larger $ -B' $) induce larger discounts on the face value, which correspond to  higher interest rates.  \n",
    "- Lower income also causes more discounting, as foreign creditors anticipate greater likelihood of default.  \n",
    "\n",
    "\n",
    "The next figure plots value functions and replicates the right hand panel of Figure 4 of [[Are08]](https://lectures.quantecon.org/zreferences.html#arellano2008default).\n",
    "\n",
    "<img src=\"https://s3-ap-southeast-2.amazonaws.com/lectures.quantecon.org/py/_static/lecture_specific/arellano/arellano_value_funcs.png\" style=\"\">\n",
    "\n",
    "  \n",
    "We can use the results of the computation to study the default probability\n",
    "$ \\delta(B', y) $ defined in [(4)](#equation-delta).\n",
    "\n",
    "The next plot shows these default probabilities over $ (B', y) $ as a heat\n",
    "map.\n",
    "\n",
    "<img src=\"https://s3-ap-southeast-2.amazonaws.com/lectures.quantecon.org/py/_static/lecture_specific/arellano/arellano_default_probs.png\" style=\"\">\n",
    "\n",
    "  \n",
    "As anticipated, the probability that the government chooses to default in the\n",
    "following period increases with indebtedness and falls with income.\n",
    "\n",
    "Next let&rsquo;s run a time series simulation of $ \\{y_t\\} $, $ \\{B_t\\} $ and $ q(B_{t+1}, y_t) $.\n",
    "\n",
    "The grey vertical bars correspond to periods when the economy is excluded from financial markets because of a past default.\n",
    "\n",
    "<img src=\"https://s3-ap-southeast-2.amazonaws.com/lectures.quantecon.org/py/_static/lecture_specific/arellano/arellano_time_series.png\" style=\"\">\n",
    "\n",
    "  \n",
    "One notable feature of the simulated data is the nonlinear response of interest rates.\n",
    "\n",
    "Periods of relative stability are followed by sharp spikes in the discount rate on government debt."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "\n",
    "<a id='arellano-ex1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "To the extent that you can, replicate the figures shown above\n",
    "\n",
    "- Use the parameter values listed as defaults in the `__init__` method of the `Arellano_Economy`.  \n",
    "- The time series will of course vary depending on the shock draws.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions\n",
    "\n",
    "Compute the value function, policy and equilibrium prices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ae = Arellano_Economy(&beta;=.953,        # Time discount factor\n",
    "                      &gamma;=2.,          # Risk aversion\n",
    "                      r=0.017,       # International interest rate\n",
    "                      &rho;=.945,        # Persistence in output\n",
    "                      &eta;=0.025,       # st dev of output shock\n",
    "                      &theta;=0.282,       # Prob of regaining access\n",
    "                      ny=21,         # Number of points in y grid\n",
    "                      nB=251,        # Number of points in B grid\n",
    "                      tol=1e-8,      # Error tolerance in iteration\n",
    "                      maxit=10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the bond price schedule as seen in figure 3 of Arellano (2008)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Create \"Y High\" and \"Y Low\" values as 5% devs from mean\n",
    "high, low = np.mean(ae.ygrid) * 1.05, np.mean(ae.ygrid) * .95\n",
    "iy_high, iy_low = (np.searchsorted(ae.ygrid, x) for x in (high, low))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 6.5))\n",
    "ax.set_title(\"Bond price schedule $q(y, B')$\")\n",
    "\n",
    "# Extract a suitable plot grid\n",
    "x = []\n",
    "q_low = []\n",
    "q_high = []\n",
    "for i in range(ae.nB):\n",
    "    b = ae.Bgrid[i]\n",
    "    if -0.35 <= b <= 0:  # To match fig 3 of Arellano\n",
    "        x.append(b)\n",
    "        q_low.append(ae.Q[iy_low, i])\n",
    "        q_high.append(ae.Q[iy_high, i])\n",
    "ax.plot(x, q_high, label=\"$y_H$\", lw=2, alpha=0.7)\n",
    "ax.plot(x, q_low, label=\"$y_L$\", lw=2, alpha=0.7)\n",
    "ax.set_xlabel(\"$B'$\")\n",
    "ax.legend(loc='upper left', frameon=False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Draw a plot of the value functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Create \"Y High\" and \"Y Low\" values as 5% devs from mean\n",
    "high, low = np.mean(ae.ygrid) * 1.05, np.mean(ae.ygrid) * .95\n",
    "iy_high, iy_low = (np.searchsorted(ae.ygrid, x) for x in (high, low))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 6.5))\n",
    "ax.set_title(\"Value Functions\")\n",
    "ax.plot(ae.Bgrid, ae.V[iy_high], label=\"$y_H$\", lw=2, alpha=0.7)\n",
    "ax.plot(ae.Bgrid, ae.V[iy_low], label=\"$y_L$\", lw=2, alpha=0.7)\n",
    "ax.legend(loc='upper left')\n",
    "ax.set(xlabel=\"$B$\", ylabel=\"$V(y, B)$\")\n",
    "ax.set_xlim(ae.Bgrid.min(), ae.Bgrid.max())\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Draw a heat map for default probability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "xx, yy = ae.Bgrid, ae.ygrid\n",
    "zz = ae.default_prob\n",
    "\n",
    "# Create figure\n",
    "fig, ax = plt.subplots(figsize=(10, 6.5))\n",
    "hm = ax.pcolormesh(xx, yy, zz)\n",
    "cax = fig.add_axes([.92, .1, .02, .8])\n",
    "fig.colorbar(hm, cax=cax)\n",
    "ax.axis([xx.min(), 0.05, yy.min(), yy.max()])\n",
    "ax.set(xlabel=\"$B'$\", ylabel=\"$y$\", title=\"Probability of Default\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot a time series of major variables simulated from the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "T = 250\n",
    "y_vec, B_vec, q_vec, default_vec = ae.simulate(T)\n",
    "\n",
    "# Pick up default start and end dates\n",
    "start_end_pairs = []\n",
    "i = 0\n",
    "while i < len(default_vec):\n",
    "    if default_vec[i] == 0:\n",
    "        i += 1\n",
    "    else:\n",
    "        # If we get to here we're in default\n",
    "        start_default = i\n",
    "        while i < len(default_vec) and default_vec[i] == 1:\n",
    "            i += 1\n",
    "        end_default = i - 1\n",
    "        start_end_pairs.append((start_default, end_default))\n",
    "\n",
    "plot_series = y_vec, B_vec, q_vec\n",
    "titles = 'output', 'foreign assets', 'bond price'\n",
    "\n",
    "fig, axes = plt.subplots(len(plot_series), 1, figsize=(10, 12))\n",
    "fig.subplots_adjust(hspace=0.3)\n",
    "\n",
    "for ax, series, title in zip(axes, plot_series, titles):\n",
    "    # Determine suitable y limits\n",
    "    s_max, s_min = max(series), min(series)\n",
    "    s_range = s_max - s_min\n",
    "    y_max = s_max + s_range * 0.1\n",
    "    y_min = s_min - s_range * 0.1\n",
    "    ax.set_ylim(y_min, y_max)\n",
    "    for pair in start_end_pairs:\n",
    "        ax.fill_between(pair, (y_min, y_min), (y_max, y_max),\n",
    "                        color='k', alpha=0.3)\n",
    "    ax.grid()\n",
    "    ax.plot(range(T), series, lw=2, alpha=0.7)\n",
    "    ax.set(title=title, xlabel=\"time\")\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "filename": "arellano.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Default Risk and Income Fluctuations"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}