insurance_incentives.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Insurance and Incentives\n",
    "\n",
    "*By Sebastian Graves and Thomas Sargent*\n",
    "\n",
    "This notebook computes optimal contracts for the three examples that lead off chapter 21 of \n",
    "**Recursive Macroeconomic Theory, Fourth edition** by Lars Ljungqvist and Thomas Sargent.\n",
    "\n",
    "The examples illustrate different sorts of tradeoffs between insurance and incentives that emerge under different\n",
    "limits on enforcement and information.  \n",
    "\n",
    "In each of the three economies, a planner or money-lender designs an efficient contract to supply insurance to a risk-averse consumer who receives an exogenous  random stream of a non-storable endowment.\n",
    "\n",
    "The only way that the consumer to smooth consumption across states and time is to interact with the planner. \n",
    "\n",
    "The three models differ in the constraints that they impose  on the planner.  \n",
    "\n",
    "These constraints express the planner's  limited ability either to enforce a contract or to observe the consumer's endowment\n",
    "\n",
    "Each of the examples uses a version of what we have nicknamed **dynamic programming squared**\n",
    "\n",
    "In a dynamic programming squared problem, a value function from one Bellman equation is an argument of another Bellman equation.  \n",
    "\n",
    "In the examples below, a planner or money lender's value function will have as an argument the value of a villager\n",
    "that satisfies a Bellman equation\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Three models of a villager and a money lender\n",
    "\n",
    "Imagine a  village with a large number of ex ante\n",
    "identical households. Each household has preferences over\n",
    "consumption streams that are ordered by\n",
    "$$ E_{-1}\\sum_{t=0}^\\infty \\beta^t u(c_t),  $$\n",
    "where $u(c)$ is an increasing, strictly concave, and twice\n",
    "continuously differentiable  function,\n",
    " $\\beta \\in (0,1)$ is a discount factor, and $E_{-1}$ is the mathematical expectation\n",
    " not conditioning on any information available at time $0$ or later.  \n",
    " \n",
    "Each household\n",
    "receives a stochastic endowment stream $\\{y_t\\}_{t=0}^\\infty$,\n",
    "where for each $t \\geq 0$, $y_t$ is independently and\n",
    "identically distributed according to the discrete\n",
    "probability distribution ${\\rm Prob} (y_t = \\overline y_s) = \\Pi_s,$\n",
    "where $s \\in \\{1, 2, \\ldots ,S\\}\\equiv {\\bf S}$ and\n",
    "$\\overline y_{s+1}>\\overline y_s$.  \n",
    "\n",
    "The consumption\n",
    "good is not storable.  \n",
    "\n",
    "At time $t \\geq 1$, the\n",
    "household has received a history of endowments\n",
    "$h_t = (y_t, y_{t-1}, \\ldots, y_0).$\n",
    "\n",
    "Endowment processes are distributed independently and identically\n",
    " both across time and\n",
    "across households.\n",
    "\n",
    "\n",
    "##### Competitive equilibrium\n",
    "\n",
    "In  this setting, if there were a competitive equilibrium with\n",
    "complete markets, at date\n",
    "$0$ households would trade history- and date-contingent claims. \n",
    "\n",
    "Since  households are  ex ante\n",
    "identical, each household would consume the per capita\n",
    "endowment in every period, and its lifetime utility would be\n",
    "\n",
    "$$ v_{\\rm pool} = \\sum_{t=0}^\\infty\n",
    "\\beta^t \\, u\\!\\left(\\sum_{s=1}^S \\Pi_s \\overline y_s\\right) =\n",
    " {1 \\over 1-\\beta}\\, u\\!\\left(\\sum_{s=1}^S \\Pi_s \\overline y_s\\right) .\n",
    "                                                         $$\n",
    "                                                         \n",
    " Households would thus insure away all\n",
    "risks from their individual endowment processes.\n",
    "\n",
    "But the\n",
    " incentive constraints that we are about to specify make\n",
    "this allocation unattainable. \n",
    "\n",
    "For each specification of incentive\n",
    "constraints, we shall solve a planning problem for an efficient\n",
    "allocation that respects those constraints.\n",
    "\n",
    "\n",
    "   Following a tradition started by\n",
    "Edward Green (1987) [*Lending and the Smoothing of Uninsurable\n",
    " Income*, in Edward C. Prescott and Neil Wallace, editors, **Contractual Arrangements for\n",
    " Intertemporal Trade**, Minnesota Studies in Macroeconomics series, Vol.\n",
    " 1, Minneapolis: University of Minnesota Press, pp. 3--25], we assume that a *moneylender* or *planner* is\n",
    "the only person in the  village who has access to\n",
    "a risk-free loan market outside the village.\n",
    "\n",
    "The moneylender can borrow or lend at a constant one-period\n",
    "risk-free gross interest rate  $R=\\beta^{-1}$.\n",
    "\n",
    "Households cannot borrow or lend with each other,\n",
    "and can  trade only with the moneylender.  \n",
    "\n",
    "Furthermore,\n",
    "we assume that the moneylender is committed   to honor  his\n",
    "promises.\n",
    "\n",
    "We will study three distinct environments in which there are three alternative types  of incentive constraints.\n",
    "\n",
    "\n",
    "**Enviroment a.** Both the money lender and the household observe   the household's  history of endowments at each  time $t$.\n",
    "Although the moneylender can commit to honor a\n",
    "contract, households  cannot commit and at any time are\n",
    "  free to walk away from an arrangement\n",
    "with the moneylender\n",
    "and live in perpetual autarky thereafter.  They must be induced not to do so\n",
    "by the structure of\n",
    "the contract.\n",
    "This is a model of *one-sided commitment* in which the\n",
    "contract must be *self-enforcing*.  That is, it must be structured to induce   the household to  prefer to\n",
    "conform to it.\n",
    "\n",
    "**Environment b.** Households *can* make commitments  and enter\n",
    "into enduring and   binding contracts with the moneylender,\n",
    "but they have private\n",
    "information about their own incomes. The moneylender\n",
    "can see neither their income nor their consumption. Instead,\n",
    " exchanges between the moneylender and a household must\n",
    "be based on the household's own reports about income\n",
    "realizations. An incentive-compatible contract  induces\n",
    "a household to report its income truthfully.\n",
    "\n",
    "**Environment c.** The environment is the same as  b except  that now households have access to a storage technology that\n",
    "cannot be observed by the moneylender.\n",
    "Households can store nonnegative amounts of goods at a risk-free\n",
    "gross return of $R$ equal to the interest rate that\n",
    "the moneylender faces in the outside credit market.\n",
    "Since the moneylender can both borrow and lend at the interest\n",
    "rate $R$ outside of the village,\n",
    "the private storage technology does not change the economy's\n",
    "aggregate resource constraint, but it does  affect the set of\n",
    "incentive-compatible contracts between the moneylender and the\n",
    "households.\n",
    "\n",
    "\n",
    "#### Preview\n",
    "\n",
    "\n",
    "When we compute efficient allocations for each of these three\n",
    "environments, we  find that the dynamics of the implied\n",
    "consumption allocations differ dramatically. \n",
    "\n",
    " \n",
    " We shall see that the dynamics\n",
    "of consumption outcomes evidently differ substantially across the\n",
    "three environments, increasing monotonically and then flattening out in environment a,\n",
    "stochastically heading south in environment b, and stochastically heading north in\n",
    "environment c.\n",
    "These sample path properties will reflect how the optimal  contracts cope with the three different frictions that we have put into the environment. \n",
    "\n",
    "Chapter 21 of RMT4  explains why  sample paths of consumption differ\n",
    "so much across these three settings."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Three computed contracts\n",
    "\n",
    "\n",
    "For all three environments discussed, consumers have a utility function:\n",
    "\n",
    "$$u(c) = - \\gamma^{-1} \\exp(-\\gamma c)$$\n",
    "\n",
    "We set $\\gamma = 0.7$, and the discount factor, $\\beta$ to 0.8. \n",
    "\n",
    "The consumers receive an iid endowment  that  can take any integer in the range $[\\bar y_1,...,\\bar y_{5}] = [6,...,10]$. \n",
    "\n",
    "The probability of each realisation is $\\Pi_s = \\frac{1-\\lambda}{1-\\lambda^{5}}\\lambda^{s-1}$ with $\\lambda = 0.4$.\n",
    "\n",
    "As mentioned above, an interesting  benchmark case is a complete markets environment.\n",
    "\n",
    "Because  all households are *ex ante* identical,  in a complete markets economy each household would consume the per capita endowment in every period, and its lifetime utility would be:\n",
    "\n",
    "$$ v_{pool} = \\frac{1}{1-\\beta} u \\left( \\sum_{s=1}^S \\Pi_s \\bar y_s \\right) =  \\frac{u(c_{pool})}{1-\\beta} $$\n",
    "\n",
    "Later we will compare the consumption paths for each environment to that which would occur in the complete markets environment. \n",
    "\n",
    "In each environment, we compute allocations  for the situation in which  the planner or money lender just breaks even.\n",
    "\n",
    "## Environment a\n",
    "\n",
    "The first environment is one in which the planner is able to commit, but households are not.\n",
    "\n",
    "At any time households are free to walk away from an arrangement with the planner, and live in perpetual autarky thereafter.\n",
    "\n",
    "RMT4 shows how this problem can be written in a recursive form. \n",
    "\n",
    "Equations 21.3.4 to 21.3.8 in RMT4 express  the planners's problem as:\n",
    "\n",
    "\\begin{align}\n",
    "\t&P(v) = \\max_{c_s,w_s} \\sum_{s=1}^S \\Pi_s \\left[ (\\bar y_s - c_s) + \\beta P(w_s) \\right] \\\\\n",
    "\t&\\text{subject to} \\\\\n",
    "\t&\\sum_{s=1}^S \\Pi_s \\left[ u(c_s) + \\beta w_s \\right] \\geq v \\\\\n",
    "\t&u(c_s) + \\beta w_s \\geq u(\\bar y_s) + \\beta v_{aut} \\text{ , s = 1,...,S} \\\\\n",
    "\t&c_s \\in [c_{min},c_{max}] \\\\\n",
    "\t&w_s \\in [v_{aut},\\bar v]\n",
    "\\end{align}\n",
    "\n",
    "where $w_s$ is the promised value with which the consumer will enter the next period, given that $y = \\bar y_s$ this period. \n",
    "\n",
    "The first constraint is a promise keeping constraint, while the second set of constraints are participation constraints. $[c_{min},c_{max}]$ is a bounded set, while $\\bar v$ just needs to be a very large number. \n",
    "\n",
    "\n",
    "The value of autarky to the households is:\n",
    "\n",
    "$$ v_{aut} = \\frac{1}{1-\\beta} \\sum_{s=1}^S \\Pi_s u(\\bar y_s) $$\n",
    "\n",
    "Below we solve the moneylender's problem in this environment by approximating $P(v)$ using Chebyshev polynomials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import minimize, fsolve\n",
    "from scipy.interpolate import UnivariateSpline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy.polynomial.chebyshev as cheb\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameter values\n",
    "gamma = 0.7\n",
    "beta = 0.8\n",
    "lamb = 0.4\n",
    "S = 5\n",
    "y_grid = np.linspace(6,5+S,S)\n",
    "prob_grid = np.zeros(S)\n",
    "for i in range(S):\n",
    "    prob_grid[i] = (1 - lamb)/(1-lamb**S)*lamb**(i)\n",
    "    \n",
    "# Utility function\n",
    "u = lambda c: -gamma**(-1)*np.exp(-gamma*c)\n",
    "u_inv = lambda u: np.log(-gamma*u)/(-gamma)\n",
    "\n",
    "# Calculate complete markets consumption    \n",
    "c_pool = np.dot(prob_grid,y_grid)\n",
    "\n",
    "# Calculate value of autarky\n",
    "v_aut = 1/(1-beta)*np.dot(prob_grid, u(y_grid))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Functions used in each environment\n",
    "\n",
    "# Nodes and basis matrix for Chebyshev approximation\n",
    "\n",
    "def Cheb_basis(order,lb,ub):\n",
    "    # Calculate roots of Chebyshev polynomial\n",
    "    k = np.linspace(order, 1, order)\n",
    "    roots = np.cos((2*k - 1)*np.pi/(2*order))\n",
    "    # Scale to approximation space\n",
    "    s = lb + (roots - -1)/2*(ub-lb)\n",
    "    # Create basis matrix\n",
    "    Phi = cheb.chebvander(roots, order-1)\n",
    "    return s, Phi\n",
    "\n",
    "# Value Function Iteration\n",
    "\n",
    "def Bellman_Iterations(s, Phi, P_fun, x_store, coeff, tolc=1e-6, bnds=None, cons=(), max_iters=100):\n",
    "    global x, c\n",
    "    c = coeff\n",
    "    order = Phi.shape[1]\n",
    "    iters = 0\n",
    "    diff = 1\n",
    "    \n",
    "    while diff > tolc:\n",
    "        # 1. Maximization, given value function guess\n",
    "        P_iter = np.zeros(order)\n",
    "        for i in range(order):\n",
    "            x = s[i]\n",
    "            res = minimize(P_fun, x_store[i], method = 'SLSQP', bounds = bnds, constraints=cons, tol=1e-15)\n",
    "            x_store[i] = res.x\n",
    "            P_iter[i] = -P_fun(res.x)\n",
    "        # 2. Bellman updating of Value Function coefficients\n",
    "        c1 = np.linalg.solve(Phi, P_iter)\n",
    "        # 3. Compute distance and update\n",
    "        diff = max(abs(c1 - c))\n",
    "        print(diff)\n",
    "        c = np.copy(c1)\n",
    "        iters = iters + 1\n",
    "        \n",
    "        if iters >= max_iters:\n",
    "            print('Convergence failed after {} iterations'.format(iters))\n",
    "            break\n",
    "            \n",
    "    if diff < tolc:\n",
    "        print('Convergence achieved after {} iterations'.format(iters))\n",
    "\n",
    "    return c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.02046457448057115\n",
      "0.01420417271952934\n",
      "0.007919316517300268\n",
      "0.004362488108637486\n",
      "0.0024521591210679983\n",
      "0.0014181717553074513\n",
      "0.0008450151371991454\n",
      "0.0005180673399127755\n",
      "0.0003258997519962614\n",
      "0.0002096164890845742\n",
      "0.00013734153609978872\n",
      "9.134953200051754e-05\n",
      "6.149551030254496e-05\n",
      "4.179812741511579e-05\n",
      "2.8630491596959295e-05\n",
      "1.973504389185532e-05\n",
      "1.367544312058655e-05\n",
      "9.5197594243146e-06\n",
      "6.6539162798529006e-06\n",
      "4.66802753196216e-06\n",
      "3.286358093546049e-06\n",
      "2.323138710569328e-06\n",
      "1.666813717604576e-06\n",
      "1.1999303339838008e-06\n",
      "8.667152331387484e-07\n",
      "Convergence achieved after 25 iterations\n"
     ]
    }
   ],
   "source": [
    "# Value Function Approximation\n",
    "# Set bounds and approximation order\n",
    "v_min = v_aut\n",
    "v_max = -0.065\n",
    "c_min = 0\n",
    "c_max = 50\n",
    "order = 70\n",
    "\n",
    "# Calculate nodes and basis matrix\n",
    "s, Phi = Cheb_basis(order, v_min, v_max)\n",
    "\n",
    "# Bounds for Maximisation\n",
    "lb = np.concatenate([np.ones(S)*c_min, np.ones(S)*v_min], axis=0)\n",
    "ub = np.concatenate([np.ones(S)*c_max, np.ones(S)*v_max], axis=0)\n",
    "\n",
    "# Initialize Value Function coefficients and goess for c,w\n",
    "y = (c_pool - u_inv(s*(1-beta)))/(1-beta)\n",
    "c = np.linalg.solve(Phi, y)\n",
    "x_init = np.concatenate([np.ones(S)*c_min, np.ones(S)*v_min], axis=0)\n",
    "\n",
    "# Function to minimize and constraints\n",
    "def P_fun(x):\n",
    "    scale = -1 + 2*(x[S:2*S] - v_min)/(v_max - v_min)\n",
    "    P = np.dot(cheb.chebvander(scale,order-1),c)\n",
    "    P_fun = - prob_grid.dot((y_grid - x[0:S]) + beta*P)\n",
    "    return P_fun\n",
    "\n",
    "def cons12(y):\n",
    "    global x\n",
    "    return prob_grid.dot(u(y[0:S]) + beta*y[S:2*S]) - x\n",
    "\n",
    "cons1 = ({'type': 'ineq', 'fun': lambda y: u(y[0:S]) + beta*y[S:2*S] - u(y_grid) - beta*v_aut},\n",
    "        {'type': 'ineq', 'fun': cons12})\n",
    "\n",
    "bnds1 = np.concatenate([lb.reshape(2*S, 1), ub.reshape(2*S, 1)], axis = 1)\n",
    "\n",
    "# Bellman Iterations\n",
    "NBell = 5\n",
    "tolc = 1e-6\n",
    "diff = 1\n",
    "iters = 1\n",
    "x_store = {}\n",
    "for i in range(order):\n",
    "    x_store[i] = x_init\n",
    "    \n",
    "c = Bellman_Iterations(s, Phi, P_fun, x_store, c, bnds=bnds1, cons=cons1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcFNW5//HP07OwKFwMIEaMQYgYWQR10LivQdyXIBJN3GK83FyXGKOX5BeuxFcWo0YlGsM1KiYGI0HBJXpdEkGjUWRQ3MB4gRhEBQcSFAeY6eX5/VHVPU0zDMV01wzT832/4NVdp05VnZrq6qfPOVWnzN0RERHZmkR7F0BERDoGBQwREYlEAUNERCJRwBARkUgUMEREJBIFDBERiUQBQzo1M5tqZpPauxxtxcx+b2anRcz7spkNjbtM0nGY7sOQjsDM3gX6Aem85Hvc/ZL2KdH2xczuAVa4+w9ayLMPcD8w1COc+GY2DjjL3b9SsoJKh1bZ3gUQ2QYnu/uf2mpjZlbp7qm22l4b+HdgepRgEXoEmGpmn3X3D2Msl3QQapKSDs/Mzjez583sRjP7l5n93cyOD+eNN7PagvxXmNkj4ft7zOxH4fsjzWyFmf2Xma0EpoXp3zSzJWb2TzN7xMx2zVuXm9kEM/u/cNu/NDPLK9cLZnazma01s2VmdnCY/p6ZfWRm5+Wtq0u4D8vNbFXYXNatoGxXhst9aGYXhPMuBs4BrjazT83s0S38qY4Hns3b3iAze8bM1pjZajObbma9svPdfSOwABjdykMjZUYBQ8rFgcDfgD7A9cBd4Rf3I8BeZrZnXt6zgfu2sJ5dgM8AnwcuNrOjgZ8C44DPAv8gaNbJdxIwChgR5juuoFyvA73Dbd4f5v0C8DXgNjPbMcz7M2AwMDKc3x/474Ky/VuY/g3gl2a2k7vfAUwHrnf3Hd395MKdMrMdgD3Cv1EuOdy3XYG9gc8BkwsWXRzul4gChnQoD4W/1LP/v5k37x/u/mt3TwO/Ifhy7+fu64GHga8ChIHjiwSBpDkZ4Bp3b3D3DQS/3O9291fcvQH4HnCQmQ3IW+Y6d1/r7suBOQRf+Fl/d/dpYblmEHwpXxuu/ymgEfhCGNy+CVzh7v9093XAT4DxeetKhssm3f1x4FNgr4h/u2zNYV02wd2XuPvTYVnqgJuAIwqWW5e3rHRy6sOQjuS0FvowVmbfuPv6sFUo+8v9PuDnwLUEtYuHwkDSnLqwKSZrV+CVvHV/amZrCH7lv1u4bWB93nYBVuW93xCuozBtR6Av0B1YEJYdghpARV7eNQV9KoXbasna8LUHsBHAzHYGfgEcFqYngH8VLNcjb1np5FTDkM7gKaCPmY0kqGlsqTkKoLBD+AOC5ikg17TTG3i/xGVcTRA8hrp7r/D/v7l71IDQYke2u9cDSwmavLJ+Gi63j7v3JGgis4JF9wZei1gGKXMKGFL2wl/lDwA3EPRPPL0Ni98HXGBmI82sC0Ez0Tx3f7fEZcwAvwZuDn/5Y2b9zey4lpfMWQUM3Eqex9m0yakHQbPWWjPrD1yVnznc3/3Ztr+XlDEFDOlIHg2vAsr+n70Ny94HHAvM3JZLZd39z8Ak4EHgQ2AQm/YrlNJ/AUuAl8zsE+BPRO+juAsYEvbtPLSFPHcA51hTm9cPgf2Aj4HHgFkF+U8B5rr7B9uwD1LGdOOeSCdiZvcBf3D3LQWV/LzzgG+4+5vxl0w6AgUMERGJRE1SIiISiQKGiIhEooAhIiKRlNWNe3369PEBAwa0dzFERDqMBQsWrHb3vlHyllXAGDBgALW1tVvPKCIiAJjZP6LmVZOUiIhEooAhIiKRKGCIiEgkZdWH0ZxkMsmKFSvYuHHj1jNLWejatSu77bYbVVVV7V0UkbJS9gFjxYoV9OjRgwEDBpA3bLSUKXdnzZo1rFixgj322KO9iyNSVsq+SWrjxo307t1bwaKTMDN69+6tGqVIDMo+YAAKFp2MjrdIPMq+SUpEtp274x48XSmTex++5r3PuJPxpvwQLOPu4WvBe5rJX5i+heWy286ECza/XFO5N1tH3vtMfv6w0E37tGle8tPJX39zf5fN15nxTctKC2WicH8Lpps7NhmH7tUVTDhiUEyfhiYKGCJl4rcvvssv5yzJTbtDOuOkMk4m45t9+QdfXs1/WUvHsnOPLgoY5aKiooLhw4fnpsePH8/EiRO3eT1Tp06le/funHvuuaUsXkmtXbuW++67j29961vNzt+wYQNjxozhmWeeoaKiotk8jY2NHHvssTzzzDNUVuojGtVfl6yhIZVhzNBdcmmVFUaFGYmEkTDDADNImIGBYSQsSDMseC3Il31vZpvmAyoSlsufyGsJbMobLJzIWy77PvzXtHyi+TLkby9/vWYtlz1/OTbZl6blKFw3YT7L7kd+ubdcjsJtZpdtWibYJs3tz2Z/35bX2dzfv63obGwD3bp1Y+HChUWvZ8KECc2mp1Kp7eaLde3atdx+++1bDBh33303Z5xxxhaDBUB1dTXHHHMMM2bM4JxzzomrqGWnvjHFHn124Lqv7NPeRZEytX18y7SRHz76Fos++KSk6xyya0+uOXloq5YdMGAA5513Ho8++ijJZJKZM2cyePBgBg4cyMKFC+nVqxcAX/jCF3jhhRf41a9+xY477sh3v/tdjjzySA4++GBeeOEFTjnlFMaOHcuFF15IXV0dffv2Zdq0aey+++6cf/759OzZk9raWlauXMn111/P2LFjmTt3Ltdccw39+vVj4cKFnHHGGQwfPpwpU6awYcMGHnroIQYNGkRdXR0TJkxg+fLlANxyyy0ccsghTJ48meXLl7Ns2TKWL1/Ot7/9bS677DImTpzI0qVLGTlyJF/+8pe54YYbNtnn6dOnc9999wHw6aefcuqpp/Kvf/2LZDLJj370I0499VQATjvtNL73ve8pYGyD+oYUO3bpVKe0tLFOcZVUe9uwYQMjR47M/Z8xY0ZuXp8+fXjllVf4j//4D2688UYSiQSnnnoqs2cHj6ueN28eAwYMoF+/fputd+3atTz77LNceeWVXHLJJZx77rm8/vrrnHPOOVx22WW5fB9++CHPP/88f/zjHzdpCnvttdeYMmUKb7zxBvfeey/vvPMOL7/8MhdddBG33norAJdffjlXXHEF8+fP58EHH+Siiy7KLf/222/z5JNP8vLLL/PDH/6QZDLJddddx6BBg1i4cOFmwaKxsZFly5aRHVG4a9euzJ49m1deeYU5c+Zw5ZVXkn0C5LBhw5g/f36Rf/nOZX1jmu7VW665iRSrU/0caW1NoFgtNUmdccYZAOy///7MmjULgLPOOotrr72WCy64gPvvv5+zzjqr2WXz01988cXc8l//+te5+uqrc/NOO+00EokEQ4YMYdWqVbn0UaNG8dnPfhaAQYMGMXr0aACGDx/OnDlzAPjTn/7EokWLcst88sknrFu3DoATTzyRLl260KVLF3beeedN1t2c1atX52pNEHSyfv/73+e5554jkUjw/vvvs2rVKnbZZRcqKiqorq5m3bp19OjRo8X1SuDThhQ7VHeqU1ramD5d7axLly5A0DGeSqUAOOigg1iyZAl1dXU89NBD/OAHP2h22R122GGL683vCMtuAyD/Ge756YlEIjedSCRyZclkMrz44ot069Zti2UvLP+WdOvWbZMb6qZPn05dXR0LFiygqqqKAQMGbDK/oaGBrl27trhOabK+Mc0OapKSGKlJajtkZpx++ul85zvfYe+996Z3795bXebggw/m/vvvB4Iv4kMPPbQkZRk9ejS33XZbbnprnfc9evTI1UAK7bTTTqTT6VxQ+Pjjj9l5552pqqpizpw5/OMfTcPyr1mzhr59+2o8qG1Q35Ciexc1SUl8FDDaQGEfRpRLas866yx+97vfbbE5qtAvfvELpk2bxj777MO9997LlClTii12br21tbXss88+DBkyhKlTp7aYv3fv3hxyyCEMGzaMq666arP5o0eP5vnnnwfgnHPOoba2lpqaGqZPn84Xv/jFXL45c+ZwwgknlGQfOoNUOkNDKqMmKYmVeRndpVNTU+OFT9xbvHgxe++9dzuVSAq9+uqr3HTTTdx7770t5jvjjDP46U9/yl577dWq7XS24/7xhiQjfvgUk04awjcO1aCLEp2ZLXD3mih5VcOQNrXvvvty1FFHkU6nt5insbGR0047rdXBojNa3xj0H+2gq6QkRqq/Spu78MILW5xfXV29Xd/Nvj2qbwgCRnd1ekuMVMMQKQP1DUGNbUd1ekuMFDBEykB92CTVXZ3eEiMFDJEykK1h6CopiZMCRhtYuXIl48ePZ9CgQQwZMoQTTjiBd955J7btzZ07l5NOOqnFPAsXLuTxxx/fpvW+++67mBmTJk3Kpa1evZqqqiouueSSbVrXPffcs03LzJ07l7/+9a/btI3OJNfprSYpiZECRszcndNPP50jjzySpUuXsmjRIn7yk59sdRiNuLUmYAAMHDiQP/7xj7npmTNnMnTotg25srU7wpujgNGyXA1Dnd4SIwWMmM2ZM4eqqqpNhiYfOXIkhx12GO7OVVddxbBhwxg+fHhuUMK5c+dyxBFHMG7cOAYPHszEiROZPn06BxxwAMOHD2fp0qUAnH/++UyYMIHDDjuMwYMHb/JFnlVfX8+FF17IqFGj2HfffXn44YdpbGzkv//7v5kxY0ZuMMTm8jWnW7du7L333mTvd5kxYwbjxo3LzX/00Uc58MAD2XfffTn22GNzgXHy5MlcfPHFjB49erMroB577DEOOuggVq9eTV1dHV/5ylcYNWoUo0aN4oUXXuDdd99l6tSp3HzzzYwcOZK//OUvzJw5k2HDhjFixAgOP/zwIo5QecheJaWAIXGK9dNlZr2AO4FhBE8UvNDdXyzIcyRwC1AFrHb3I8L0K4CLwuXeAC5w940U438nwso3ilrFZnYZDsdft8XZb775Jvvvv3+z82bNmsXChQt57bXXWL16NaNGjcp9+b322mssXryYz3zmMwwcOJCLLrqIl19+mSlTpnDrrbdyyy23AEEz0bPPPsvSpUs56qijWLJkySbb+PGPf8zRRx/N3Xffzdq1aznggAM49thjufbaa6mtrc0N+/H973+/2XzNjVc1fvx47r///twggbvuuisffPABAIceeigvvfQSZsadd97J9ddfz89//nMAFixYwPPPP0+3bt245557AJg9ezY33XQTjz/+ODvttBNnn302V1xxBYceeijLly/nuOOOY/HixUyYMCE3tDsEAyQ++eST9O/fn7Vr10Y9WmUr2+ndrUpNUhKfuH+OTAGecPexZlYNdM+fGQaU24Ex7r7czHYO0/sDlwFD3H2Dmf0BGA/cE3N529Tzzz/PV7/6VSoqKujXrx9HHHEE8+fPp2fPnpFGkgUYN24ciUSCPffck4EDB/L2229vso2nnnqKRx55hBtvvBGAjRs35p5tESVfc3dLjxkzhkmTJtGvX7/Nhi5ZsWIFZ511Fh9++CGNjY3ssUfTXcennHLKJoMYzpkzh9raWp566il69uwJtDw6br5DDjmE888/n3HjxuVG/O3M6htSdKuqoCLRdk9fk84ntoBhZj2Bw4HzAdy9EWgsyHY2MMvdl4d5PiooWzczSxIEmg+KLlQLNYG4DB06lAceeKDZeS0NyxJlJFlgs8czFk67Ow8++OBmd03PmzcvUr7mVFdXs//++/Pzn/+ct956i0cffTQ379JLL+U73/kOp5xyCnPnzmXy5Mm5eYW1lYEDB7Js2TLeeecdamqCkQlaGh0339SpU5k3bx6PPfYYI0eOZOHChZEGaSxX9RqpVtpAnH0YA4E6YJqZvWpmd5pZYfvGYGAnM5trZgvM7FwAd38fuBFYDnwIfOzuTzW3ETO72Mxqzay2rq4uvr1ppaOPPpqGhgZ+/etf59Lmz5/Ps88+y+GHH86MGTNIp9PU1dXx3HPPccABB2zT+mfOnEkmk2Hp0qUsW7Zssy/84447jltvvTUXnF599VVg81Flt5RvS6688kp+9rOfbfYl/fHHH9O/f38AfvOb37S4js9//vPMmjWLc889l7feegvY8ui4heVdunQpBx54INdeey19+vThvffea3Fb5W59Q0pXSEns4gwYlcB+wK/cfV+gHigcprUS2B84ETgOmGRmg81sJ+BUYA9gV2AHM/tacxtx9zvcvcbda/r27RvTrrSemTF79myefvppBg0axNChQ5k8eTK77rorp59+Ovvssw8jRozg6KOP5vrrr2eXXXbZpvXvtddeHHHEERx//PFMnTp1s+dHTJo0iWQyyT777MOwYcNyl8QeddRRLFq0KNfpvaV8WzJ06FDOO++8zdInT57MmWeeyWGHHUafPn0ilX/69OmceeaZLF26dIuj45588snMnj071+l91VVXMXz4cIYNG8bhhx/OiBEjov7JytKnDWndtCexi220WjPbBXjJ3QeE04cBE939xLw8E4Gu7j45nL4LeCKcPcbdvxGmnwt8yd2/1dI2O9toteeffz4nnXQSY8eObe+ibHfK+bg35+xfv0QynWHmhIPbuyjSwWwXo9W6+0rgPTPLtpEcAywqyPYwcJiZVZpZd+BAYDFBU9SXzKy7BY3yx4TpItKM+kbVMCR+cX/CLgWmh1dILQMuMLMJAO4+1d0Xm9kTwOtABrjT3d8EMLMHgFeAFPAqcEfMZe1wspemitQ3pOjfS4+zlXjFGjDcfSFQWNWZWpDnBuCGZpa9BrimROXY7OohKV/l9FCwqNY3pDSOlMSu7O/07tq1K2vWrOmUXyKdkbuzZs2azTr/y50uq5W2UPafsN12240VK1awPV5yK/Ho2rUru+22W3sXo824O/UNKbrraXsSs7IPGFVVVZvcbSxSbhrTGVIZVw1DYlf2TVIi5W597lkYqmFIvBQwRDq4T/U8b2kjChgiHdz6xuzzvBUwJF4KGCIdXK6GoSYpiZkChkgH1/R4VtUwJF4KGCIdXO7xrLpxT2KmgCHSwTU9nlVNUhIvBQyRDi7bJKXBByVuChgiHVy9rpKSNqKAIdLB1TekMIOuVTqdJV76hIl0cPUNaXaortSIzBI7BQyRDm59o57nLW1DAUOkg/tUz8KQNqKAIdLBrW9M0101DGkDChgiHVy9ahjSRhQwRDq4+saUhgWRNqGAIdLBrW9Ia+BBaRMKGCId3KcNKd20J21CAUOkg1vfmNawINImFDBEOjB3p74xxY66SkragAKGSAe2IZnGXY9nlbahgCHSgTU9C0M1DImfAoZIB6an7UlbUsAQ6cCanuetgCHxU8AQ6cDWh8/C0OCD0hYUMEQ6sKbHs6qGIfFTwBDpwJo6vRUwJH4KGCIdWH3ued5qkpL4KWCIdGDrwyYpDQ0ibUEBQ6QDqw87vfU8DGkLChgiHVh9Q4rKhFFdoVNZ4hfrp8zMepnZA2b2tpktNrODmslzpJktNLO3zOzZbVlWpLOrbwiehWFm7V0U6QTibvicAjzh7mPNrBronj/TzHoBtwNj3H25me0cdVkRCZqkNCyItJXYAoaZ9QQOB84HcPdGoLEg29nALHdfHub5aBuWFSlrmYyTyjjpjJPKZMJXb3pNO3XrGjTwoLSZOD9pA4E6YJqZjQAWAJe7e31ensFAlZnNBXoAU9z9txGXBcDMLgYuBth9991j3B0pB39fXc/yf64nlc6QyjipdPBlnEw7qXSGZCZ4TWecZNpJZ+dlMmHepnzptJPMfpGnm77UkwXrTGectAd50nlf+Nky5AeCZF5+92j7NGrATvH+0URCcQaMSmA/4FJ3n2dmU4CJwKSCPPsDxwDdgBfN7KWIywLg7ncAdwDU1NREPMWkszrl1udZF16Kui0qE0ZFwqhMGJUVCaoqstNN76sqEkFaRYKqMH/36srccomEhXkTVBhUViRy680tm7edRF56hRmVFXnpiQSJMN+Iz/WK4S8lsrk4A8YKYIW7zwunHyD40i/MszqsOdSb2XPACOAvEZYV2SbpjLOuIcX4UZ/jqwfsnvuirqwIvoArKyz3viI/LWHqVBYhxoDh7ivN7D0z28vd/0ZQi1hUkO1h4DYzqwSqgQOBmyMuK7JNkukMALv37q5f5SKtEHdv2aXA9PAqp2XABWY2AcDdp7r7YjN7AngdyAB3uvubW1o25rJKmcsGDN2zINI6sQYMd18I1BQkTy3IcwNwQ8RlRVotlQ66uKoUMERaJdKZY2ZnmNn/mdnHZvaJma0zs0/iLpxIKWVrGAoYIq0TtYZxPXCyuy+OszAicWrMBQx1YIu0RtSfWqsULKSjS4ZNUtWVqmGItEbUGkatmc0AHgIasonuPiuWUonEQE1SIsWJGjB6AuuB0XlpDihgSIfRmFLAEClGpIDh7rqkVTq8pPowRIoS9Sqp3cxstpl9ZGarzOxBM9st7sKJlFJSl9WKFCXqmTMNeATYFegPPBqmiXQY6sMQKU7UM6evu09z91T4/x6gb4zlEik5XVYrUpyoAWO1mX3NzCrC/18D1sRZMJFSS6rTW6QoUc+cC4FxwErgQ2BsmCbSYeg+DJHiRL1KajlwSsxlEYmV+jBEitNiwDCzq939ejO7leC+i024+2WxlUykxNSHIVKcrdUwssOB1MZdEJG4aXhzkeK0GDDc/dHw7Xp3n5k/z8zOjK1UIjFQp7dIcaKeOd+LmCay3crduKdOb5FW2VofxvHACUB/M/tF3qyeQCrOgomUmvowRIqztT6MDwj6L04BFuSlrwOuiKtQInHIXSWVUA1DpDW21ofxGvCamd0HGPBFgqul/ubujW1QPpGSSaWdyoSRSKiGIdIaUYc3/zLwP8BSgsCxh5n9u7v/b2wlEymxZDqjDm+RIkQNGDcBR7n7EgAzGwQ8BihgSIfRmM6o/0KkCFF/bn2UDRahZcBHMZRHJDbJdEbDgogUIWoN4y0zexz4A0EfxpnAfDM7A/SoVukYkimnUh3eIq0WNWB0BVYBR4TTdcBngJPRo1qlg0imM1RVqklKpLX0iFbpNBrV6S1SlEgBw8z2AC4FBuQv4+4awVY6jGQ6o3GkRIoQtUnqIeAugkezZuIrjkh8kmlXDUOkCFEDxkZ3/8XWs4lsv5K6rFakKFEDxhQzuwZ4CmjIJrr7K7GUSiQGjSn1YYgUI2rAGA58HTiapiYpD6dFOoRkOsMOXaJ+5EWkUNSz53RgoMaPko5MfRgixYl69rwG9IqzICJxUx+GSHGi1jD6AW+b2Xw27cPQZbXSYeg+DJHiRA0Y18RaCpE2oPswRIoT9U7vZ1uzcjPrBdwJDCPoJL/Q3V8syHMkcAtQBax29yPy5lUQPMDpfXc/qTVlEMlKptSHIVKMqHd6ryP4wgeoJvhyr3f3nltZdArwhLuPNbNqoHvBensBtwNj3H25me1csPzlwGKCR8KKFEVjSYkUJ9LPLXfv4e49w/9dga8At7W0jJn1BA4nuEMcd29097UF2c4GZrn78jDPR3nL7wacSFBDESma+jBEitOqs8fdH2Lr92AMJBjVdpqZvWpmd5rZDgV5BgM7mdlcM1tgZufmzbsFuJqtDEViZhebWa2Z1dbV1W3jnkhnkkq7+jBEihC1SeqMvMkEUENTE1VL694PuNTd55nZFGAiMKkgz/7AMUA34EUze4kgkHzk7gvCPo4tcvc7gDsAampqtlYm6cSS6QyVuqxWpNWiXiV1ct77FPAucOpWllkBrHD3eeH0AwQBozDPanevB+rN7DlgBEGgOcXMTiB4FkdPM/udu38tYnlFNpHJOKmMOr1FihHb8zDcfaWZvWdme7n73whqEYsKsj0M3GZmlQSd6QcCN7v7TOB7kLuK6rsKFlKMZCZo2VTAEGm9SGePmV1vZj3NrMrM/mxmq80syhf4pcB0M3sdGAn8xMwmmNkEAHdfDDwBvA68DNzp7m+2bldEtiyZDlor1Ych0npRm6RGu/vVZnY6QTPSmcAc4HctLeTuCwn6O/JNLchzA3BDC+uYC8yNWE6RZiVT2RqG+jBEWivqz62q8PUE4Pfu/s+YyiMSi2Q6DBiVqmGItFbUGsajZvY2sAH4lpn1BTbGVyyR0mpMqw9DpFhRb9ybCBwE1Lh7Eqhn61dJiWw31IchUrxteZrM3sCA8IqmrN+WuDwisUiqhiFStKg37t0LDAIWAukw2VHAkA6iUZ3eIkWLWsOoAYa4u+6klg5Jnd4ixYt69rwJ7BJnQUTipD4MkeJFrWH0ARaZ2cvoiXvSAakPQ6R4UQPG5DgLIRK3pstq1Ych0lqRn7hnZv2AUWHSy/nPrhDZ3jXd6a0ahkhrRR1LahzBWE9nAuOAeWY2Ns6CiZRSrg9Dnd4irRa1Ser/AaOytYrwTu8/EQxZLrLdy/ZhVCbUJCXSWlF/biUKmqDWbMOyIu1OQ4OIFC9qDeMJM3sS+H04fRbweDxFEim9lJqkRIrWYsAwsy8A/dz9qvAxrYcCBrwITG+D8omUhC6rFSne1s6eW4B1AO4+y92/4+5XENQubom7cCKlktRltSJF21rAGODurxcmunstMCCWEonEQH0YIsXb2tnTtYV53UpZEJE4JVNBH4YChkjrbe3smW9m3yxMNLNvAAviKZJI6SXTGSoSRoUuqxVpta1dJfVtYLaZnUNTgKgBqoHT4yyYSCkl0xn1X4gUqcWA4e6rgIPN7ChgWJj8mLs/E3vJREqoMZ1Rc5RIkaKOJTUHmBNzWURik0xnNLS5SJF0BkmnkEy5ahgiRdIZJJ1CMp2hqlJ9GCLFUMCQTkF9GCLF0xkknYL6MESKpzNIOoVkWn0YIsXSGSSdQjKdoVL3YYgURQFDOoXGlPowRIqlM0g6BfVhiBRPZ5B0CkEfhpqkRIqhgCGdQlKX1YoUTWeQdArBjXv6uIsUQ2eQdArJtKsPQ6RIsZ5BZtbLzB4ws7fNbLGZHdRMniPNbKGZvWVmz4ZpnzOzOeEyb5nZ5XGWU8qfhjcXKV6k0WqLMAV4wt3Hmlk10D1/ppn1Am4Hxrj7cjPbOZyVAq5091fMrAewwMyedvdFMZdXypT6MESKF9sZZGY9gcOBuwDcvdHd1xZkOxuY5e7RUa/HAAALyklEQVTLwzwfha8fuvsr4ft1wGKgf1xllfKn+zBEihfnGTQQqAOmmdmrZnanme1QkGcwsJOZzTWzBWZ2buFKzGwAsC8wr7mNmNnFZlZrZrV1dXWl3QMpG8m0U61Ob5GixHkGVQL7Ab9y932BemBiM3n2B04EjgMmmdng7Ewz2xF4EPi2u3/S3Ebc/Q53r3H3mr59+8awG1IO1IchUrw4A8YKYIW7Z2sGDxAEkMI8T7h7vbuvBp4DRgCYWRVBsJju7rNiLKeUuUzGSWU0+KBIsWI7g9x9JfCeme0VJh0DFHZaPwwcZmaVZtYdOBBYbGZG0Pex2N1viquM0jkkMxkABQyRIsV9ldSlwPTwCqllwAVmNgHA3ae6+2IzewJ4HcgAd7r7m2Z2KPB14A0zWxiu6/vu/njM5ZUylEw7gO7DEClSrAHD3RcCNQXJUwvy3ADcUJD2PKAGZymJZCpbw9BHSqQY+sklZS+ZDgJGpWoYIkXRGSRlrzEMGGqSEimOziApe9k+jKpKNUmJFEMBQ8petklKV0mJFCfuq6Q6hv+dCCvfaO9SSEz6N6a4v/pjBv+lByyobu/iiJTeLsPh+Oti34x+cknZy3jQJKUGKZHiqIYBbRKZpf0s/vs/Gf8/L/K70Qdy6J592rs4Ih2WahhS9lJp3YchUgoKGFL2spfV6hGtIsXRGSRlT0ODiJSGziApe7qsVqQ0dAZJ2UuqD0OkJBQwpOw1plTDECkFnUFS9nJ9GOr0FimKziApe+rDECkNnUFS9tSHIVIaChhS9hpVwxApCZ1BUvaSqXB4cwUMkaLoDJKyl0xnSBhUJNQkJVIMBQwpe8l0RrULkRLQWSRlrzGd0bAgIiWgs0jKXjKd0cCDIiWgs0jKXjLluqRWpAQUMKTsqQ9DpDR0FknZUx+GSGnoLJKypxqGSGnoLJKyl0o7VZXqwxAplgKGlL1G1TBESkJnkZQ9NUmJlIbOIil7ybSr01ukBHQWSdkLahjqwxAplgKGlL3GlJqkREpBZ5GUPfVhiJSGziIpe8m0hgYRKYVYA4aZ9TKzB8zsbTNbbGYHNZPnSDNbaGZvmdmzeeljzOxvZrbEzCbGWU4pb6phiJRGZczrnwI84e5jzawa6J4/08x6AbcDY9x9uZntHKZXAL8EvgysAOab2SPuvijm8koZ0mi1IqUR21lkZj2Bw4G7ANy90d3XFmQ7G5jl7svDPB+F6QcAS9x9mbs3AvcDp8ZVVilvjSmNJSVSCnHWMAYCdcA0MxsBLAAud/f6vDyDgSozmwv0AKa4+2+B/sB7eflWAAc2txEzuxi4GGD33XdvVUGfXrQKA6orE1RVJKiutPA1nK4IXqsqjOrKIL26IoGZ2sU7AvVhiJRGnAGjEtgPuNTd55nZFGAiMKkgz/7AMUA34EUzewlo7uz25jbi7ncAdwDU1NQ0m2drLv39K2xMZrZ5uaoKCwNJUxDpUpk3XRlMZ183CT55QalLXv7qigTVlRWbLds0L9G0zXC6a1WCLpUVdKlMkNBzqzejPgyR0ogzYKwAVrj7vHD6AYKAUZhndVjrqDez54ARYfrn8vLtBnwQV0Ef+s9DSKacxnSGxlSGZPjamA7eJ9O+WXpz+RpSm6c3pjJ82pCiIZkhmWmal0w7yVSGhnC6VCoT1hRgcgGnIghmVWFQq2qa7lKQ3iUMUF0qK4L5uWBVsUkA2yRvwXJVFbbd1L4yGSeVcQUMkRKILWC4+0oze8/M9nL3vxHUIgo7rR8GbjOzSqCaoNnpZuBtYE8z2wN4HxhP0N8Riy/u0jOuVUfi7kFQygaa7P90mo3JzQNRMp0hmfFcvoZUNl86nM5sMi+73oZUhoZkho83JJvm5dKDfA2pDN6qeloTM3K1rWww6bqFYJStGeUHt65VQf6u4XK5Zas2zZ99H+Rp2kZ+sEpmgmBcrU5vkaLFfZXUpcD08AqpZcAFZjYBwN2nuvtiM3sCeB3IAHe6+5sAZnYJ8CRQAdzt7m/FXNZ2Y2ZUVwb9I3Rp37LkB6+GZDoXfBoKA0ze+43J/PSmALQxG4SSGTZuEpwy/Ku+cZO82XwNYc2sGE3BKpHr7K5UU51I0cyL/Tm5HampqfHa2tr2LoYUKZ3xXK1pYxh4NiYzeTWlTQPXxmSGDY1pNiQLglgYhDLuXHL0FxjUd8f23jWR7Y6ZLXD3mih5465hiGyzioTRvbqS7tXtXRIRyaeGXRERiUQBQ0REIlHAEBGRSBQwREQkEgUMERGJRAFDREQiUcAQEZFIFDBERCSSsrrT28zqgH+0cvE+wOoSFqcj6Iz7DJ1zvzvjPkPn3O9t3efPu3vfKBnLKmAUw8xqo94eXy464z5D59zvzrjP0Dn3O859VpOUiIhEooAhIiKRKGA0uaO9C9AOOuM+Q+fc7864z9A59zu2fVYfhoiIRKIahoiIRKKAISIikXT6gGFmY8zsb2a2xMwmtnd54mJmnzOzOWa22MzeMrPLw/TPmNnTZvZ/4etO7V3WUjOzCjN71cz+GE7vYWbzwn2eET5CuKyYWS8ze8DM3g6P+UHlfqzN7Irws/2mmf3ezLqW47E2s7vN7CMzezMvrdlja4FfhN9vr5vZfsVsu1MHDDOrAH4JHA8MAb5qZkPat1SxSQFXuvvewJeA/wz3dSLwZ3ffE/hzOF1uLgcW503/DLg53Od/Ad9ol1LFawrwhLt/ERhBsP9le6zNrD9wGVDj7sOACmA85Xms7wHGFKRt6dgeD+wZ/r8Y+FUxG+7UAQM4AFji7svcvRG4Hzi1ncsUC3f/0N1fCd+vI/gC6U+wv78Js/0GOK19ShgPM9sNOBG4M5w24GjggTBLOe5zT+Bw4C4Ad29097WU+bEmeOR0NzOrBLoDH1KGx9rdnwP+WZC8pWN7KvBbD7wE9DKzz7Z22509YPQH3subXhGmlTUzGwDsC8wD+rn7hxAEFWDn9itZLG4BrgYy4XRvYK27p8LpcjzmA4E6YFrYFHenme1AGR9rd38fuBFYThAoPgYWUP7HOmtLx7ak33GdPWBYM2llfZ2xme0IPAh8290/ae/yxMnMTgI+cvcF+cnNZC23Y14J7Af8yt33Beopo+an5oRt9qcCewC7AjsQNMcUKrdjvTUl/bx39oCxAvhc3vRuwAftVJbYmVkVQbCY7u6zwuRV2Spq+PpRe5UvBocAp5jZuwTNjUcT1Dh6hc0WUJ7HfAWwwt3nhdMPEASQcj7WxwJ/d/c6d08Cs4CDKf9jnbWlY1vS77jOHjDmA3uGV1JUE3SSPdLOZYpF2HZ/F7DY3W/Km/UIcF74/jzg4bYuW1zc/Xvuvpu7DyA4ts+4+znAHGBsmK2s9hnA3VcC75nZXmHSMcAiyvhYEzRFfcnMuoef9ew+l/WxzrOlY/sIcG54tdSXgI+zTVet0env9DazEwh+dVYAd7v7j9u5SLEws0OBvwBv0NSe/32Cfow/ALsTnHRnunthh1qHZ2ZHAt9195PMbCBBjeMzwKvA19y9oT3LV2pmNpKgo78aWAZcQPADsWyPtZn9EDiL4IrAV4GLCNrry+pYm9nvgSMJhjFfBVwDPEQzxzYMnrcRXFW1HrjA3Wtbve3OHjBERCSazt4kJSIiESlgiIhIJAoYIiISiQKGiIhEooAhIiKRVG49i4gUMrPeBIO8AewCpAmG4wBY7+4Ht0vBRGKky2pFimRmk4FP3f3G9i6LSJzUJCVSYmb2afh6pJk9a2Z/MLN3zOw6MzvHzF42szfMbFCYr6+ZPWhm88P/h7TvHog0TwFDJF4jCJ7HMRz4OjDY3Q8guAv70jDPFIJnNowCvhLOE9nuqA9DJF7zs2P3mNlS4Kkw/Q3gqPD9scCQYBQHAHqaWY/wuSUi2w0FDJF45Y9blMmbztB0/iWAg9x9Q1sWTGRbqUlKpP09BVySnQgHDhTZ7ihgiLS/y4AaM3vdzBYBE9q7QCLN0WW1IiISiWoYIiISiQKGiIhEooAhIiKRKGCIiEgkChgiIhKJAoaIiESigCEiIpH8f9wlW8Cjr5zQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Time Series Simulation\n",
    "T = 100\n",
    "np.random.seed(2)\n",
    "y_series = np.random.choice(y_grid,size = T,p = prob_grid)\n",
    "c_series = np.zeros(T)\n",
    "w_series = np.zeros(T)\n",
    "resid_series = np.zeros(T)\n",
    "pval_series = np.zeros(T)\n",
    "\n",
    "# Initialize v such that P(v) = 0\n",
    "v_find = lambda v: cheb.chebvander(-1 + 2*(v - v_min)/(v_max - v_min),order-1).dot(c)\n",
    "x = fsolve(v_find,v_max)\n",
    "\n",
    "res = minimize(P_fun,x_init,method = 'SLSQP',bounds = bnds1,constraints = cons1,tol=1e-15)\n",
    "c_series[0] = res.x[np.where(y_grid == y_series[0])[0][0]]\n",
    "w_series[0] = res.x[S + np.where(y_grid == y_series[0])[0][0]]\n",
    "\n",
    "# Simulate\n",
    "for t in range(1,T):\n",
    "    x = w_series[t-1]\n",
    "    res = minimize(P_fun, x_init,method = 'SLSQP',bounds = bnds1, constraints = cons1, tol=1e-15)\n",
    "    c_series[t] = res.x[np.where(y_grid == y_series[t])[0][0]]\n",
    "    w_series[t] = res.x[S + np.where(y_grid == y_series[t])[0][0]]\n",
    "    \n",
    "plt.plot(c_series, label = 'Environment (a)')\n",
    "plt.plot(np.ones(T)*c_pool, label = 'Complete Markets')\n",
    "plt.ylabel('Consumption')\n",
    "plt.xlabel('Time')\n",
    "plt.legend(loc = 'best');\n",
    "plt.title('Environment (a)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above simulation is equivalent to Figure 21.2.1.a in RMT. \n",
    "\n",
    "The discussion in RMT4 confirms that the household's consumption ratchets upwards over time. \n",
    "\n",
    "The consumption level is constant after the first time that the household receives the highest possible endowment.\n",
    "\n",
    "## Environment b\n",
    "\n",
    "The second environment is one in which households *can* make commitments to enter into binding contracts with the planner, but they have private information about their incomes. \n",
    "\n",
    "Consequently, incentive compatability constraints are required to ensure that households truthfully report their incomes. \n",
    "\n",
    "Equations 21.5.1 to 21.5.5 in RMT4 express the planners's problem. \n",
    "\n",
    "\\begin{align}\n",
    "\t&P(v) = \\max_{b_s,w_s} \\sum_{s=1}^S \\Pi_s \\left[ -b_s + \\beta P(w_s) \\right] \\\\\n",
    "\t&\\text{s.t.} \\\\\n",
    "\t&\\sum_{s=1}^S \\Pi_s \\left[ u(\\bar y_s + b_s) + \\beta w_s \\right] = v \\\\\n",
    "\t& C_{s,k} \\equiv u(\\bar y_s + b_s) + \\beta w_s - [ u(\\bar y_s + b_k) + \\beta w_k ] \\geq 0 \\hspace{2mm} \\forall \\hspace{2mm} s,k \\in S \\times S\\\\\n",
    "\t&b_s \\in [a - \\bar y_s,\\infty ] \\\\\n",
    "\t&w_s \\in [- \\infty, v_{max}]\n",
    "\\end{align}\n",
    "\n",
    "Here $b_s$ is the transfer that the moneylender gives to a household who reports income $y_s$ if their promised value was $v$. \n",
    "\n",
    "The promise keeping constraint remains, while the participation constraint has been replaced by a large set of incentive compatibility constraints. \n",
    "\n",
    "RMT4 shows that we can discard many of the incentive compatibility constraints.\n",
    "\n",
    "In solving the model below, we keep only the local upward and downward incentive compatibility constraints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0356136025755589\n",
      "0.030371315540371313\n",
      "0.02407399182513359\n",
      "0.01913266149864512\n",
      "0.015222370598273471\n",
      "0.012114987856719495\n",
      "0.009631147461220735\n",
      "0.007672223008313495\n",
      "0.006105497493145151\n",
      "0.00486238930683669\n",
      "0.003874043097987112\n",
      "0.0030852738234585786\n",
      "0.002458861972733928\n",
      "0.00195938541418883\n",
      "0.0015613577112887356\n",
      "0.0012446135038572947\n",
      "0.0009919573701395734\n",
      "0.0007906335097516148\n",
      "0.0006302442740278025\n",
      "0.0005023133549428849\n",
      "0.0004003561312018178\n",
      "0.0003190925702796221\n",
      "0.0002542816361241762\n",
      "0.0002026263867875855\n",
      "0.0001614498820146082\n",
      "0.00012861680815490217\n",
      "0.0001024510164100434\n",
      "8.159488184134034e-05\n",
      "6.497036918773347e-05\n",
      "5.17247087685746e-05\n",
      "4.116982454860363e-05\n",
      "3.2760247570706724e-05\n",
      "2.6062380570124333e-05\n",
      "2.072756330306902e-05\n",
      "1.6479524738599594e-05\n",
      "1.3097959921992697e-05\n",
      "1.0406276047092433e-05\n",
      "8.26447011803566e-06\n",
      "6.560744409966901e-06\n",
      "5.205745019054575e-06\n",
      "4.1285380945055294e-06\n",
      "3.27246526410363e-06\n",
      "2.5923429660679176e-06\n",
      "2.052268719410222e-06\n",
      "1.6235840192280193e-06\n",
      "1.283467490509338e-06\n",
      "1.0137738115645334e-06\n",
      "8.000368652005818e-07\n",
      "Convergence achieved after 48 iterations\n"
     ]
    }
   ],
   "source": [
    "# Set bounds and approximation order\n",
    "b_min = -20\n",
    "b_max = 20\n",
    "w_min = -150;\n",
    "w_max = -0.04;  \n",
    "v_min = -150;\n",
    "v_max = -0.04;\n",
    "v_pool = u(c_pool)/(1-beta)\n",
    "order = 70\n",
    "\n",
    "# Calculate nodes and basis matrix\n",
    "s, Phi = Cheb_basis(order,v_min,v_max)\n",
    "\n",
    "# Bounds for Maximisation\n",
    "lb = np.concatenate([np.ones(S)*b_min,np.ones(S)*w_min], axis=0)\n",
    "ub = np.concatenate([np.ones(S)*b_max,np.ones(S)*w_max], axis=0)\n",
    "\n",
    "# For initial guess, use upper bound given in RMT:\n",
    "cbar = np.zeros(order)\n",
    "upper = np.zeros(order)\n",
    "for i in range(order):\n",
    "    cbar[i] = u_inv((1-beta)*s[i])\n",
    "    upper[i] = np.dot(prob_grid,(y_grid - cbar[i])/(1-beta))\n",
    "c = np.linalg.solve(Phi,upper)\n",
    "\n",
    "# Function to minimize and constraints\n",
    "def P_fun2(x):\n",
    "    scale = -1 + 2*(x[S:2*S] - v_min)/(v_max - v_min)\n",
    "    P = np.dot(cheb.chebvander(scale,order-1),c)\n",
    "    P_fun = - prob_grid.dot(-x[0:S] + beta*P)\n",
    "    return P_fun\n",
    "\n",
    "def cons23(y):\n",
    "    global x\n",
    "    return prob_grid.dot(u(y_grid + y[0:S]) + beta*y[S:2*S]) - x\n",
    "\n",
    "cons2 = ({'type': 'ineq', 'fun': lambda x:  u(y_grid[1:S] + x[1:S]) + beta*x[S+1:2*S] - u(y_grid[1:S] + x[0:S-1]) - beta*x[S:2*S-1]},\n",
    "         {'type': 'ineq', 'fun': lambda x:  u(y_grid[0:S-1] + x[0:S-1]) + beta*x[S:2*S-1] - u(y_grid[0:S-1] + x[1:S]) - beta*x[S+1:2*S]},\n",
    "         {'type': 'eq',  'fun': cons23})\n",
    "\n",
    "bnds2 = np.concatenate([lb.reshape(2*S,1),ub.reshape(2*S,1)], axis = 1)\n",
    "\n",
    "x_store = {}\n",
    "for i in range(order):\n",
    "    x_store[i] = np.concatenate([np.zeros(S),np.ones(S)*s[i]], axis=0)\n",
    "\n",
    "c = Bellman_Iterations(s, Phi, P_fun2, x_store, c, tolc, bnds=bnds2, cons=cons2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Time Series Simulation\n",
    "T = 800\n",
    "np.random.seed(2)\n",
    "y_series = np.random.choice(y_grid,size = T+1, p = prob_grid)\n",
    "c_series = np.zeros(T)\n",
    "w_series = np.zeros(T)\n",
    "\n",
    "# Initialize v such that P(v) = 0\n",
    "v_find = lambda v: cheb.chebvander(-1 + 2*(v - v_min)/(v_max - v_min),order-1).dot(c)\n",
    "x = fsolve(v_find,v_aut)\n",
    "\n",
    "x_init = np.concatenate([np.zeros(S),np.ones(S)*x],axis=0)\n",
    "res = minimize(P_fun2,x_init,method = 'SLSQP',bounds = bnds2, constraints = cons2,tol=1e-10)\n",
    "c_series[0] = y_series[0] + res.x[np.where(y_grid == y_series[0])[0][0]]\n",
    "w_series[0] = res.x[S  + np.where(y_grid == y_series[0])[0][0]]\n",
    "x_init = res.x\n",
    "\n",
    "# Simulate\n",
    "for t in range(1,T):\n",
    "    x = w_series[t-1]\n",
    "    res = minimize(P_fun2,x_init,method = 'SLSQP',bounds = bnds2,constraints = cons2,tol=1e-10)\n",
    "    c_series[t] = y_series[t] + res.x[np.where(y_grid == y_series[t])[0][0]]\n",
    "    w_series[t] = res.x[S + np.where(y_grid == y_series[t])[0][0]]\n",
    "    x_init = res.x\n",
    "\n",
    "# Plot    \n",
    "plt.plot(c_series, label = 'Environment (b)')\n",
    "plt.plot(np.ones(T)*c_pool, label = 'Complete Markets')\n",
    "plt.ylabel('Consumption')\n",
    "plt.xlabel('Time')\n",
    "plt.title('Environment (b)')\n",
    "plt.legend(loc = 'best');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This simulation reported in the graph above confirms that in environment **b** the incentive compatibility constraints induce  the planner to introduce a downward tilt into consumption paths.\n",
    "\n",
    "## Environment c\n",
    "\n",
    "The third environment is the same as in (b), except for the additional assumption that households have access to a storage technology. \n",
    "\n",
    "A household  can store nonnegative amounts that cannot be observed by the planner.\n",
    "\n",
    "The text of RMT4 chaper 21 shows that the solution to this problem is the same as in an economy in which  each household can lend *or borrow* at the risk-free gross interest rate R, subject to the natural debt limit.\n",
    "\n",
    "Thus, the planner enables the household to relax the no-borrowing constraint implied by  the restriction that it can store only nonnegative amounts \n",
    "\n",
    "We can find the natural debt limit by iterating forward on the households budget constraint:\n",
    "\n",
    "\\begin{equation}\n",
    "\tc + k' = y + Rk\n",
    "\\end{equation}\n",
    "This iteration gives:\n",
    "\\begin{equation}\n",
    "\tk = \\frac{1}{R} \\sum_{j=0}^\\infty \\frac{c - y}{R^j}\n",
    "\\end{equation}\n",
    "\n",
    "Imposing non-negativity on consumption:\n",
    "\n",
    "\\begin{equation}\n",
    "\tk \\geq  - \\frac{1}{R} \\sum_{j=0}^\\infty \\frac{y}{R^j}\n",
    "\\end{equation}\n",
    "\n",
    "Finally, the natural debt limit is found by choosing the lowest possible value of the endowment, so that for any possible endowment stream the household can always pay back its debts:\n",
    "\n",
    "\\begin{equation}\n",
    "\tk \\geq  - \\frac{1}{R} \\sum_{j=0}^\\infty \\frac{\\bar y_{min}}{R^j} = - \\frac{\\bar y_{min}}{R-1} \\equiv \\phi\n",
    "\\end{equation}\n",
    "\n",
    "A recursive presentation of the household's problem is then:\n",
    "\\begin{align}\n",
    "\t&V(k,y) = \\max_{c,k'} u(c) + \\beta E [V(k',y')] \\\\\n",
    "\t&\\text{s.t.} \\\\\n",
    "\t&c + k' = y + Rk \\\\\n",
    "\t& k' \\geq \\phi\n",
    "\\end{align}\n",
    "\n",
    "As income is iid, we can re-write the household's problem with only one state. \n",
    "\n",
    "Define a = k + y. \n",
    "\n",
    "Then\n",
    "\\begin{align}\n",
    "\t&V(a) = \\max_{c,k'} u(c) + \\beta E [V(Rk' + y')] \\\\\n",
    "\t&\\text{subject to} \\\\\n",
    "\t&c + k' = a \\\\\n",
    "\t& k' \\geq \\phi\n",
    "\\end{align}\n",
    "\n",
    "Below we solve this latter problem using Value Function Iteration, again with Chebyshev polynomials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Update parameter values\n",
    "# Set bounds and approximation order\n",
    "R = 1/beta\n",
    "k_min = - y_grid[0]/(R - 1)\n",
    "k_max = 100\n",
    "a_min = R*k_min + min(y_grid)\n",
    "a_max = R*k_max + max(y_grid)\n",
    "order = 150\n",
    "\n",
    "# Calculate nodes and basis matrix\n",
    "s, Phi = Cheb_basis(order,a_min,a_max)\n",
    "\n",
    "# Create bounds\n",
    "bnds3 = np.array([[k_min,k_max]])\n",
    "\n",
    "# Value function\n",
    "def P_fun3(kprime):\n",
    "    global x,c\n",
    "    # Function to minimize\n",
    "    scale = -1 + 2*(R*kprime + y_grid - a_min)/(a_max - a_min)\n",
    "    P_fun = -(u(x - kprime) + beta * prob_grid.dot(cheb.chebval(scale, c)))\n",
    "    return P_fun\n",
    "\n",
    "# Initialize guess and VF coefficients\n",
    "c = np.zeros(order)\n",
    "\n",
    "x_store = {}\n",
    "for i in range(order):\n",
    "    x_store[i] = k_min\n",
    "\n",
    "c = Bellman_Iterations(s, Phi, P_fun3, x_store, c, bnds=bnds3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Time Series Simulation\n",
    "T = 800\n",
    "np.random.seed(2)\n",
    "y_series = np.random.choice(y_grid, size = T+1, p = prob_grid)\n",
    "\n",
    "a_series = np.zeros(T+1)\n",
    "c_series = np.zeros(T)\n",
    "\n",
    "# Initialise at v_aut\n",
    "def k_find(k): \n",
    "    scale = -1 + 2 * (R * k + y_grid - a_min)/(a_max - a_min)\n",
    "    return prob_grid.dot(cheb.chebval(scale,c)) - v_aut\n",
    "k0 = fsolve(k_find,0)\n",
    "a_series[0] = k0 + y_series[0]\n",
    "\n",
    "# Simulate\n",
    "for t in range(T):\n",
    "    x = a_series[t]\n",
    "    res = minimize(P_fun3, k_min,method='SLSQP',  bounds=bnds3,tol=1e-15)\n",
    "    c_series[t] = a_series[t] - res.x\n",
    "    a_series[t+1] = R * res.x + y_series[t+1]\n",
    "\n",
    "# Plot\n",
    "plt.plot(c_series, label = 'Environment (c)')\n",
    "plt.plot(np.ones(T)*c_pool, label = 'Complete Markets')\n",
    "plt.ylabel('Consumption')\n",
    "plt.xlabel('Time')\n",
    "plt.title('Environment (c)')\n",
    "plt.legend(loc = 'best')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the introduction of a storage technology for the household means that the consumption path now has an upward trend. \n",
    "\n",
    "This occurs because our parameter values satisfy  $\\beta R = 1$."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}