assignment3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Assignment 1: International Macroeconomics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import norm\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problem 1\n",
    "4) Tauchen's Method\n",
    "\n",
    "Reference: http://www.econ.nyu.edu/user/violante/NYUTeaching/Macrotheory/Spring14/LectureNotes/lecture6_14.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "σ_ε = 0.01                        # Variance of ε\n",
    "α = 0.7                           # Autoregression parameter\n",
    "m = 3                             # Tauchen multiple\n",
    "F = norm(loc=0, scale=σ_ε).cdf    # Distribution of ε\n",
    "T = 10                            # Size of Y grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_max = m * np.sqrt(σ_ε**2 / (1 - α**2))\n",
    "y_min = -y_max\n",
    "\n",
    "y_grid = np.linspace(y_min, y_max, T) # Discrete y grid\n",
    "\n",
    "d = (y_max - y_min) / (T - 1)         # Step size\n",
    "π = np.empty((T, T))                  # Transition matrix\n",
    "\n",
    "for i in range(T):\n",
    "    π[i,     0] = F(y_grid[0] - α * y_grid[i] + d / 2)\n",
    "    π[i, T - 1] = 1 - F(y_grid[T - 1] - α * y_grid[i] - d / 2)\n",
    "    for j in range(1, T - 1):\n",
    "        z = y_grid[j] - α * y_grid[i]\n",
    "        π[i, j] = F(z + d / 2) - F(z - d / 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y_grid = np.exp(y_grid)\n",
    "β = 0.96\n",
    "ρ = -np.log(β)\n",
    "r = 0.02\n",
    "# r = 1 / np.exp(-ρ) - 1       # Such that β(1 + r) = 1\n",
    "N = 100                      # Size of A grid\n",
    "\n",
    "B_bar = -np.min(Y_grid) / r  # Natural debt limit\n",
    "A_max = 15\n",
    "A_grid = np.linspace(B_bar, A_max, N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "def iterate_V(Y_grid,   # Output grid\n",
    "              A_grid,   # Asset grid\n",
    "              π,        # Transition matrix of Y\n",
    "              r=r,      # Interest rate\n",
    "              ρ=ρ,      # Risk-free rate\n",
    "              σ=3):     # Intertemporal substitution of elasticity\n",
    "    \n",
    "    β = np.exp(-ρ)                      # Discount factor\n",
    "    u = lambda c: c**(1 - σ) / (1 - σ)  # Utility function\n",
    "    \n",
    "    T = len(Y_grid)        # Size of output grid\n",
    "    N = len(A_grid)        # Size of asset grid\n",
    "    \n",
    "    # Initial guess of value function\n",
    "    # Columns = asset choices, rows = random y shocks\n",
    "    V = np.zeros((T, N))        \n",
    "    \n",
    "    V_new = np.zeros_like(V)    # To store maximised value of V\n",
    "    A_policy = np.zeros_like(V) # To store optimal choice of A'\n",
    "    \n",
    "    max_iter = 1000\n",
    "    dist = 1e3      # Initial distance\n",
    "    tol = 1e-5      # Tolerance\n",
    "    i = 0\n",
    "    \n",
    "    while i < max_iter and dist > tol:\n",
    "        \n",
    "        for y_i, Y in enumerate(Y_grid):\n",
    "            π_temp = π[y_i, :]                      # Transition given current Y\n",
    "            \n",
    "            for a_i, A in enumerate(A_grid): \n",
    "                C = (1 + r) * A + Y - A_grid        # Consumption given A and Y\n",
    "                C[C <= 0] = 1e-3                    # Non-negativity constraint\n",
    "                V_temp = u(C) + β * π_temp @ V      # Values for all possible A' choices\n",
    "                A_index = np.argmax(V_temp)         # Maximising index of A\n",
    "                A_policy[y_i, a_i] = A_index\n",
    "                V_new[y_i, a_i] = V_temp[A_index]        # Maximised value of V\n",
    "            \n",
    "        dist = np.linalg.norm(V - V_new)\n",
    "        V = V_new.copy()\n",
    "        i += 1\n",
    "        \n",
    "    if i < max_iter:\n",
    "        print(f'Converged in {i} iterations')\n",
    "        \n",
    "    return V, A_policy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Converged in 467 iterations\n"
     ]
    }
   ],
   "source": [
    "V_star, A_policy = iterate_V(Y_grid, A_grid, π)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hUVR7G8e9vZhJ67z2UgAUUIdKCCtJRxIKKK1JEinXBAigiKFJEd11kVZogiAIiIihVqkqv0iH0QCihhTIh9ewfGRTd0JJM7tyb3+d58iRzZ5J5R8P7nJw591wxxqCUUspZXFYHUEoplfG03JVSyoG03JVSyoG03JVSyoG03JVSyoE8VgcAKFy4sAkJCbE6hlJK2cr69etPGmOKpHZfQJR7SEgI69atszqGUkrZiogcvNp9Oi2jlFIOpOWulFIOpOWulFIOdN1yF5FxInJCRLZecaygiPwsIhG+zwV8x0VEPhGRPSKyWURq+DO8Ukqp1N3IyP1LoPnfjvUBFhljQoFFvtsALYBQ30dX4POMiamUUupmXLfcjTG/AKf/drg1MMH39QTg4SuOTzQpVgH5RaRERoVVSil1Y9I6517MGHMUwPe5qO94KSDyiscd9h37PyLSVUTWici66OjoNMZQSimVmox+Q1VSOZbqnsLGmNHGmDBjTFiRIqmuwb+utQdO868Fu0hISk7T9yullFW8SckM3BtF5KV4v/z8tJb78cvTLb7PJ3zHDwNlrnhcaSAq7fGubcPBM4xYvIf4RC13pZR9/HbmPA3X7OTTQydYdOqcX54jreU+C+jg+7oDMPOK4+19q2bqADGXp2/8weNOiZ+YrBccUUoFvnOJSby+M5I2m/aSFH+JHkcmUGrRZL8813W3HxCRyUADoLCIHAb6A0OBb0WkM3AIeNz38DlAS2AP4AU6+SHzHzyulFmgJC13pVSAm38yht67DnMiPoFGxzbTJu9/yVvyFKe31PPL81233I0xT13lrkapPNYAL6Y31I1y+8o9MVmnZZRSgSk6PoG3I44w88RZQkig17ExVCu2iEvevJzZ8CCPvz7cL88bEBuHpZWO3JVSgcoYw/TjZ+gXcYQLSck8cGQ1Dxf6jBxFL3J2353UDO9PyIN3+u35bV3uf4zck7TclVKB48ileHrtOsyi0+eonOTl5dOfUaXkSmIvFOTczuY81vMDv2ewdbl73DpyV0oFjmRjmBh1ivf3RpFkDI8e/oUHi44iqGA8Z3eFcV/rf1H4odKZksXW5e526WoZpVRg2HZ6D6/tPMymuPxUTTjEU+fHE1JqM96YYiQdbc5jL7yTqXlsXe46566Uslp8UgLDts5l1OniBJkgupjPuMe1GJM7iNPb6tD06U/IV7BQpueydbnrahmllJXWHN/GqzsPsCe5LDUTtnHfhkN4L5Vkd0JXKocU4/GXO1z/h/iJrctdR+5KKSvEJsYyYPM8JsWUITd56Lh3PjkjY0nETcjFYJ56/xWyZc9uaUZbl/ufI3ctd6VU5lh6ZAOvR0Rz2JSn3qUN1N5wjNiEeIok5KJajYrUb/Oo1REBm5e7x/eGqo7clVL+FhN3nrd+/5nvL4RQSDx02TUPz9E4DB4qeINo/+EbVkf8C1uXu65zV0plhh8PrqTvPi8nqMD9F1dxx6ZTXEpKoHhcLmo1qcFdjf7vhH3L2brcdZ27UsqfomNP88bvS5kXW4ESnOOFrXNJPhmHywRTKSGYdkMDa7R+JVuXu66WUUr5y2drZzD8fAHOU44mMb9y++YYLiUnUupSLu59rAFVat1tdcRrsnW562oZpVRGizixn56bV7DOfTtlkg/RetdS8pxwkS05mBAJou0HgTtav5Kty11XyyilMkpSUhLDVn3L2LgyxLtCeTB6CT3WVSKvNGTP+a3c9kZLipcvb3XMG2brctfVMkqpjLAxcjtv7NrCVvetVEzeS/fVR2l4PowzCSfZXmwzzT942eqIN83W5a4jd6VUeiQkJjBg+RS+TqqEcZWjTdRCXt5UlexyO9svrKfWgKepVjRt13i2mq3L/c85d31DVSl1c37bt5439x0gwl2N2xJ38OKqc9T21uZkwjEOhx6m6Ys9rI6YLrYud13nrpS6Wd64WN5a8S3Tza0Eu4rR7uACum+/CzdF2XJxNQ2HPU/2XLmsjpluti53XeeulLoZc3f+Sv8jZzjkupO74jfz0opE7oyvy7G4w1y6W2jR7nWrI2YYW5e7zrkrpW5ETOw5Xl85ndncQW5JokvEfJ7bU4tEk8jvl1bS5KN/Epwtm9UxM5Sty11XyyilrufbzQsYHJ3EMddd1PGu45/LsxOaXI/DsfsJalaEBx7sZXVEv7B1uevIXSl1NcdOHuO1dXNYHFydQpyi5/bFPHWoJnFJsWxMWk6rEX2sjuhXti53XS2jlErN2CXf8XFyfk4HV+f+C+t5bU1x8sffwSHvHgq2vYVW9Z1d7GDzcteRu1LqSgejDtLz9yWsyF6dEkTxwroZ5Fy/kaWSmyLly/LMf4dZHTHT2Lrc/xi561JIpbK8jxd8w+eeklzIVo3mZ5ZTddZKJPYs2TzFaNWvJ+VuqWp1xExl63LXkbtSaueBXby6ay0bgqtSNukgHVfMJ9u2LbgkLyWq1qRtv3etjmgJW5e7iOB2ia6WUSoLSkxMZMjCrxkfXJGEoMo8HL2USj8sR5LiyB5cnCc/6E/hkmWsjmkZW5c7pIzedeSuVNaybsdGekXuZnu2OwlN3EOrJatw792JWwoQUrcWD/ewx7a8/mT7cve4RFfLKJVFxMXH8c7CSUzJfjviKcsTUQspM2s5LpLJkaMEz3z8AXkKFLQ6ZkCwfbnryF2prGHJxuW8ffIYe3PUpGr8DpouWIX78F7cUpBbm9WjWafuVkcMKLYvd4/OuSvlaBe9F+m9dDI/ZL+TbO6itNs/j2LzV+LCRa7cpeg04mOy5cxpdcyAk65yF5GewHOAAbYAnYASwBSgILABeMYYE5/OnFfldrl05K6UQ01d8B0fubIRmSOMmrG/03DOaiT6EB5XYWo80YJ7HnnS6ogBK83lLiKlgFeA24wxsSLyLdAWaAl8bIyZIiIjgc7A5xmSNhUel+g6d6UcJvrkcd5YPoMFee4mL+fouOMHiizbgkgQeQqWpePw/xAUHGx1zICW3mkZD5BDRBKAnMBR4H7gH777JwAD8GO5u11Cgr6hqpRjjJ09iRHZC3M8bx3CL6yj0fwETOK95Miem1rd61C1bn2rI9pCmsvdGHNERD4CDgGxwAJgPXDWGJPoe9hhoFRq3y8iXYGuAGXLlk1rDDxunXNXygkOHtpPr80LWZbrbgonn+D5zT9ReHsdPAnnyFP6F556b4DVEW0lPdMyBYDWQHngLDANaJHKQ1NtXmPMaGA0QFhYWJrbWVfLKGV//5k5jtG5y3EmZ03uj1lJ/XnZMdQjR9xKGrzRigq3P2p1RNtJz7RMY2C/MSYaQES+B+oB+UXE4xu9lwai0h/z6nTOXSn72hmxjd4Ra1idtwYlk47w8ro55Ntfj6C4k+SvsprH+/S1OqJtpafcDwF1RCQnKdMyjYB1wBKgDSkrZjoAM9Mb8lp0tYxS9pOQkMCwn77gy3y34s1ejZanfyVsXgGMqw7ZE37lgQHtKF7uCatj2lp65txXi8h3pCx3TAQ2kjLNMhuYIiLv+459kRFBr0bPUFXKXtb9voa3j+5iU/46lEs8QJeVv5Erqi7BCUcpWnMrD73c3+qIjpCu1TLGmP7A3/9P7ANqpefn3gydc1fKHhISEnjnxzFMzX8nCcFVeOT4EqouLI24apHdLOGRD1+kQJHiVsd0DD1DVSnld0tXL2FAzAl2FqhH5YQIHl62l2ynapEt4SBl74+haYeBVkd0HNuXu47clQpcF70X6TN7HDML3Y3LU4YnDy+i8tJKuFzVye5ZyBMjepMrbz6rYzqS7cs9yO0iNiHJ6hhKqb/5bvYM/hUM+wvfQ7W4bTRfepE8JyrhSYgk9JGi3PPoYKsjOprty11H7koFltOnTtNr2WTm5q9NTrx02DeXNh9PwOUSIppW4rHB3xEUnM3qmI5n+3LX1TJKBY7xM77hv3nycqRAOGGxm2g7eQY11+4moryHCv0/pG2d1M5zVP5g+3J3u4REPYlJKUsdORJF77WzWJTvbvITQ5eds3hixGRis8HWJ6rz6Dtf4fbYvm5sxfb/tXVvGaWsNWLaOMYULMWJfHWof2EtbcdP5fbtR9hRJYg73x9JrWr1rI6YJdm+3N0ul5a7UhbYvWc3fXcs49fCd1M0+TgvbZrGI6O/JyYX7Hz2Ph7tNdLqiFma7cvdo2+oKpXphk4dxYQioZzNVYNGMSv5x+eTqHDwNFvvyEG9DyZRr/xtVkfM8mxf7m49iUmpTLNu+a+8G7OHtUVrUyrpMJ1X/cQDk+YSnR/2vfIQj7/wgdURlY/tyz1l5K6rZZTyp7i4OAZO/YyppWvizV6NB6J/4R/DvyR7XCybw/LQ5F/TKVisjNUx1RVsX+46clfKv5bOm8tgTrG5TEPKJ+6nyZZpzCi0lOMd4flbe/Nk/Y5WR1SpsH2565y7Uv5x6aKXvj+M4vsStUmiII8cW8L2mC/4sSC0Sg6l79MTyZUzj9Ux1VXYvtzdLpderEOpDDbrh2l8lNOwu2RDqiTsJnzzAn4ssppQl4v3qg3kvpqtrY6orsP25e5x68hdqYxy/uxZes+fwI9F6uIhkbaRP7Pi0jgWFHLzhKs6vTt+QbBuHWALti93nXNXKmNMnjKB4YXycaDofdwZt5WqW2ezoMhWbic7r9f6mLDbG1gdUd0E25e7rpZRKn1OHImiz6rvmV+0Djm5SPt9c1iUNIkVBYUO2e7h1Xaf4nK7rY6pbpLty93tEpINJCcbXC6xOo5StvLJl8P5snQlogrWp5Z3AyV3zmBu4QPUSMxFnwYjubVCTasjqjSyfbl7fIWeZAwutNyVuhG7tm/h3d2/saTsPeQ3Z2i/awbrg6YRmRdeyNOKbu0G62jd5mxf7m6XC4CkZEOQ/i4qdV1Dx33EpHJ3cDJfXe49v5rn102joWsta/M1JfSZEeQvrNcxdQLbl/vlkbuumFHq2tavXc6gYztYUb4xxZKP8crW7+lz8hNOuArx+31jubvh41ZHVBnI9uXuvjwto2vdlbqq/uM/ZGq5MM7luoumZ5bzyqYJhLl2sbrIo9ze/t8Uz1vA6ogqg9m+3D3uyyN3XTGj1N8tWzyPYZeOsz6kCaWTIumyeR6vnRtJpLsk25tNpXad5lZHVH5i+3L/Y+Su0zJK/eHSpUv0//a/TC9dm0vZb+PB6GX03jKS8q5jrCzVnrvaDaFMztxWx1R+ZPty1zl3pf7qxx+nMTwoga1lGlMhcR//2PQrL10cx76gEPY/NIu61e+xOqLKBLYv9ytXyyiVlXnPn+PNWaOZWaIeybhoE7WIN7aPoITnPCvLP0/YP94lSLcOyDJsX+46clcKvpk0gs+LliCiZGNujd9JlzU/cW/iSk7lLAePfUzdW2pYHVFlMtuX+59z7vqGqsp6zkQf481Fk5ldMpwg4nn60ALm7rqF9z3t6NVsIM/UDfnj34jKWmxf7jpyV1nV6HH/YmzZShwq1pDql7Zw59aTTD11O/eEFmbwI9UoUzCn1RGVhWxf7pdHJYm6zl1lEVEHInhr3Xx+DrmPXFyg0965TD9QnShPUT5scxttapZGREfrWZ3ty/3yOnd9Q1VlBR+PGcrECndwtFB9ans3UHrrJSafuYPmtxflvYdvp2ie7FZHVAEiXeUuIvmBsUBVwADPAruAqUAIcAB4whhzJl0pr+HyahmdllFOtmn9CoZGbmZZxaYUNKfpsuMnJh+uwZGcwXz+9O20qFbC6ogqwLjS+f3DgXnGmFuAO4EdQB9gkTEmFFjku+03Hj2JSTnce6MG8fTZiyzLW4t7z6+l+qqDLL1Qj8dqlGbhq/dqsatUpXnkLiJ5gXuBjgDGmHggXkRaAw18D5sALAV6pyfktfwx566rZZTDLFs8h/94o1hZ+QGKJx+lx+759O7+JujlS9UNSM/IvQIQDYwXkY0iMlZEcgHFjDFHAXyfi6b2zSLSVUTWici66OjoNIfw6BuqyoH6jh1CV3KyOuddND39G9PKlE0pdqVuUHrm3D1ADeBlY8xqERnOTUzBGGNGA6MBwsLC0tzMureMcpIff/iGz7IlsrFiC8okHeLFvWt4pVsvq2MpG0pPuR8GDhtjVvtuf0dKuR8XkRLGmKMiUgI4kd6Q1+LRN1SVA8R6vfT7dgQzyoYTTzCtopfyXvhjlGj8kNXRlE2ludyNMcdEJFJEqhhjdgGNgO2+jw7AUN/nmRmS9Cr0DFVld5O/Hs3YwnnYVq4ZFRP30jHyEF2e7Wl1LGVz6V3n/jLwtYgEA/uATqTM438rIp2BQ4BfL+/y537uOnJX9nLhXAxvzhrDjyXDMQhtohYxtNVz5G6Sz+poygHSVe7GmE1AWCp3NUrPz70ZOueu7OiLccMZX7oUe0o15rb4HXQ5eZ6nnn7N6ljKQex/hqqullE2cvRIJO/8Np15IfUIJp52B+cz8ImXyZFT94FRGcv25a4jd2UXn4waxqSKt3CoaANqXPqd5+OCaNXRb6eAqCzO9uWuq2VUoIvYvY2BW5awMLQReThP571zGfScrllX/mX7ctfVMiqQfTByCJNDq3OsYH3qXlxHj5wluU+LXWUC25e77ueuAtH6Vb8w7PgOllVpQaHkaF7cPYd+3d6yOpbKQtK7cZjl3Lrlrwow744eTPuL8fyS524axqzk6/y5tNhVptORu1IZZPGCmQxPOMXq0JaUSI6iZ8R8enXTKRhlDduXu66WUVaL9Xp5b/Jwppevx0VPKZqf+pUBdzUhpFFLq6OpLMz25R50ebWMrnNXFpj1/SQ+y2HYVKEFZZMO8sreNbzU7Q2rYyll/3J3uQQRXS2jMles18vb0/7LjDL1SCCY1seX8H7DJynSWDdbV4HB9uUOKfPuOueuMsvXk0Yxtkg+dpRtSmjCHjpFHeHZjrrRlwosjih3t0t0zl353elTJ+k3bwI/lQxHMDxxeCGDW3chd17d6EsFHkeUu8fl0pG78qvRX3zMhLLl2FuyEVXjtvPc2Yu0feZ1q2MpdVWOKHcduSt/OXokkn7Lv2de+fpk4xLtD8zn3Sd1oy8V+BxR7ilz7vqGqspY/xn1AZMq3sbhIvdRM/Z3XkzOQctOutGXsgdHlLuO3FVG2rl9E4N2LmdRaGPyco7n9szl/S56MpKyF0eUu8clus5dZYghIwczJbQGxwuEU+/COl7NU4b6WuzKhhxR7m63jtxV+qxZsZgPT+7h1yotKZx8gpd3zaFvd90PRtmXI8pdV8uo9BgwejBTK9bmbO4w7o9ZyVshNajaSItd2Zsjyl3n3FVaLJz3A58knWFNaEtKJh3huYj5vNZdp2CUMzii3HW1jLoZsV4v7075hOkh9fBSmhanfqX/XU0IafyA1dGUyjCOKHcduasbNWPal3yex8Pm8s0pl3iAnvvX8EJXPRlJOY8jyl33llHXE+v18uZ3nzKzdDiJeHj4+BIGNnySIk0etjqaUn7hiHLXkbu6lq8mfsYXxQuxs0wTKidE0CnqGJ10oy/lcI4od4/Lpevc1f85feokby+YyE+lw3GTxJOHf2ZQ66660ZfKEhxR7jpyV3/3+ZiPmBhSkf3F76da3Da6xcTT5hm9iIbKOhxR7h63EJeYZHUMFQAiD+6j/5qfmF+xATnx0mH/PAa0fUU3+lJZjiPKXUfuCuDfI4fydaWqHCl8L3d7N/ISeWj2bB+rYyllCUeUu66Wydp2bNnA+3tWsbhyE/ITQ5eIuQzsqicjqazNZXWAjKAj96xr0MjBPHkimkX56xF+YQPjg4wWu1I4ZuSue8tkNat+XcRHZ/bxW5WWFE0+ziu7ZvNW975Wx1IqYKR75C4ibhHZKCI/+W6XF5HVIhIhIlNFJDj9Ma9NR+5ZS7/RQ+iYICzPXYNGZ1cwtWgRLXal/iYjpmX+Cey44vYHwMfGmFDgDNA5A57jmnRvmaxh/uzvaDV7PGNCW5Ar2csbu3/m60de4NZqNayOplTASVe5i0hp4AFgrO+2APcD3/keMgHw+/ndbpeQpCcxOVas10uv8UN5MUcxNuSoxgMnf+GH0Nt4tbuuhFHqatI75/4foBeQx3e7EHDWGJPou30YKJXaN4pIV6ArQNmyZdMVwuPW1TJO9d3U8YzKF8yWkOaUT9xPh4P76P7ca1bHUirgpXnkLiIPAieMMeuvPJzKQ1NtXWPMaGNMmDEmrEiRImmNAeicuxNdOBdDj4kf8kaRKuwKrsSjxxYzu8Y9WuxK3aD0jNzDgYdEpCWQHchLykg+v4h4fKP30kBU+mNem8flIiFJ59ydYvyXIxhfshi7yzThloRddD52imfav2p1LKVsJc3lbox5E3gTQEQaAK8bY54WkWlAG2AK0AGYmQE5r0lH7s4QfTyKfkumMrtsOEEk8I/IBQxq85JuHaBUGvhjnXtvYIqIvA9sBL7ww3P8hZ6han+fjf6IieUrcaBYQ+68tJXnvYaH2/eyOpZStpUh5W6MWQos9X29D6iVET/3RunI3b4O7N3Juxt/ZkGlBuTkIp32zeWdp/6po3Wl0skhZ6imjNyNMaSsxlR28NHIIXxT6Q6iCt1DLe8GegQV4v7OunWAUhnBEeXudqUs+kk24NZuD3hbN61m8IENLKncjPzmDN0j5jCg61tWx1LKURyxcZjH1+h6lmrgGzRyMG1PxbA4X13qX1jLhByixa6UHzhk5J5S7jrvHrh+Wzqff5+PZEWVlhRLPkaPXXPo011LXSl/cUS5e1yXR+5a7oGo35jBTKtYl5hcd9H4zArevqUetzTSYlfKnxxR7n+M3HV/mYAye+YUPguKY32llpROiqTbnjX07N7b6lhKZQmOKHcduQeWWK+X/lNHML1cPeLIzoPRyxgY/iglGreyOppSWYYjyv3yahmdc7felCljGJs/F1tDmlExcR8dDh2ka+eeVsdSKstxRLn/OXLX1TJWuXAuhrdmjmFWqXAMQpuoRbzXvAMFmxS2OppSWZIjyl1Xy1hr3JfDGV+yJBGlG3Nr/E6ei47h6Xa6e6NSVnJEuf+5zl3LPTNFH4/i7SVTmVM2nCDidaMvpQKII8pdR+6Z77+jPuSrCpU5WKwh1S9t4YVY4SHd6EupgOGIcv9jzl2XQvpdxO5tDNqyhJ9DG5KLC3TeO5dBz+l+MEoFGkeUu66WyRzDRg1hcqU7OVqwPnUurqdHjmI00GJXKiA5otx1bxn/2rR+BUMjN7MstBkFzWmej5hDf90PRqmA5oyNw3TO3W8GjhpMu7MXWZa3FveeX8vEXMFa7ErZgCNG7m49QzXDLVs8h/94j7KyckuKJx+lx+759O6uUzBK2YUjyt2jc+4Zqu/YIXxXoS7nc1an6enf6FetIaGNWlgdSyl1ExxR7jpyzxg//vANn2VLZGPFFpRJOsSLe9fwSjdd3qiUHTmi3P+cc9c3VNMi1uul37cjmFE2nHiCaXViKe/Vf4wSjR+yOppSKo0cUe5uXeeeZpO/Hs2YwnnYXq4ZFRP30jHyEF2e1Y2+lLI7R5T75aWQOud+4y6ci6HPj2P5sUQ4AI8fWciQh7qQu0k+i5MppTKCM8pd59xvyhfjhjOuTGn2lmzEbfE76HLyPE+1e93qWEqpDOSIctczVG/M0SORvPPbdOaF1COYeNodnM/AJ17Wjb6UciBHlLuO3K/vk1HD+KriLUQWbUCNS7/zYnw2Huiol7xTyqkcUe5uXS1zVTu3b2LwzuUsDG1EHs7rRl9KZRGOKHcduafug5FD+Cb0Lo4XCKfuxXX0yFmS+7TYlcoSHFHuup/7X61f9QvDju/gl8rNKGRO8eLuOfTrpvvBKJWVOGTjsJSXoevc4d1Rg2l/MZ5f8txNg3OrmZQ/lxa7UlmQM0buus6dxQtmMjzhFKsrt6REchQ9I+bTq5tOwSiVVTmi3LPynHus18t7k4czvXw9LnpK0fz0bwyo3piQRi2tjqaUslCay11EygATgeJAMjDaGDNcRAoCU4EQ4ADwhDHmTPqjXl1WXS3zw/SvGJkTNlVoQbnEg7yybw0vdXvD6lhKqQCQnjn3ROA1Y8ytQB3gRRG5DegDLDLGhAKLfLf9yi1Za+Qe6/Xy6sRhvFqgAtuyVaH18SX8dMfdWuxKqT+keeRujDkKHPV9fV5EdgClgNZAA9/DJgBLAb+eLeNyCS7JGnPuX038jHHFCrKjTFNCEyJ49uhROnXQjb6UUn+VIXPuIhIC3AWsBor5ih9jzFERKXqV7+kKdAUoW7ZsujN4XC5Hj9xPnzpJv/kT+bF0PVwk8+SRnxn0UFdy59WNvpRS/y/d5S4iuYHpQA9jzDnxTZFcjzFmNDAaICwsLN2t7HaJY0fuo8b+mwnlQthX4n6qxm2n69lYnminUzBKqatLV7mLSBApxf61MeZ73+HjIlLCN2ovAZxIb8gb4XGJ49a5Rx7cx7trfmJehXvIziXaH5jPu0/qRl9KqetLz2oZAb4Adhhj/n3FXbOADsBQ3+eZ6Up4g9xuIdFBq2X+M+oDJlW8jcOF7yUsdhMvJOekZSfd6EspdWPSM3IPB54BtojIJt+xt0gp9W9FpDNwCHg8fRFvjMcljphz37l9E4N2LmdhaGPyEUPXiLm811VPRlJK3Zz0rJb5DbjaBHujtP7ctHK7hCSbT8sMGTmYyaE1OVEgnPoX1tIjT1nqa7ErpdLAEWeogr1Xy6xZsZgPT+7h1yotKZJ8gpd2z+Ft3Q9GKZUOjin3lNUy9ptzHzB6MFMr1uZs7po0OruStyvV5tZGWuxKqfRxTLnbbc594bwf+CTpDGtCW1Iy6QjPRSzgte46BaOUyhiOKXe7rHOP9Xp5d8onTA+ph5fStDz1C+/e/SBlGj9gdTSllIM4qtwDfeQ+Y9qXfJ7Hw+byzQlJPEDP/Wt4oevrVsdSSjmQY8rd4w7ckXus18ub333KzNLhJOHmkWOLGdS0PQWbPGx1NKWUQ4XQuRwAAAjvSURBVDmm3N0Bulpm4sRP+aJ4EXaVaULlhN10ijpOp46vWh1LKeVwjil3T4Ctljl96iRvL5jIT6XDcZNE28ifGdLmRd06QCmVKRxT7u4A2lvm8zEfMTGkIvuL30+1uG10i4mnTXvd6EsplXkcU+4elxCfaO3IPfLgPvqv+Yn5FRuQEy8d98+jf9tXdLSulMp0jil3q1fL/HvkEL6uVI0jhe/lbu9G/ukuQONn/X4RKqWUSpVjyt1j0Tr3HVs28P6e1Syu3JT8xNB1zxze66JnmCqlrJWea6gGFCtWywwaOZgnTpxkUf661L+wnvFBRotdKRUQHDZyz5w591W/LuKjs/v4rUpLiiYf55Vds3mre99MeW6llLoRjin3lIt1+H/k3m/0EKZVqk1Mrho0OruCtyvV4dZGWuxKqcDimHL395z7/NnfMUIusC60BaWSDtNtz1p6dtcrIymlApNjyt1f69xjvV76T/2E78uFE0sID578hf61HqRM4wcz/LmUUiqjOKbc/TFy/3bKF4zJn50tIc0pn7ifDgf30f251zL0OZRSyh8cU+4ZuVrmwrkY+s4azcyS4STj4tFji3m/aXsKNnkkQ36+Ukr5m2PKPaNWy4yf8AnjS5Rgd6km3JKwi87HTvFMe93oSyllL84p93Sulok+HkW/JVOZXSacIBL4R+QCBrV5SbcOUErZknPKPR1z7p+O/pCJ5UM5WKwhd17ayvNew8Pte2VwQqWUyjyOKfe0zLkf2LuTdzf+zIJKDcnJRTrtm8s7T/1TR+tKKdtzTLnf7Mj9o5FD+KbSHUQVuofa3g38M6gQ93fWC1QrpZzBMeV++QLZxhhE5KqP27ppNYMPbGBJ5WYUMGfoHjGHAV11PxillLM4ZuMwjyul0K81en9/1GDanjrH4nx1uefCWr7MIVrsSilHcs7I3Z1S7onJBo/7r/f9tnQ+/z4fyYrKLSmWfIweu+bQp7uWulLKuRxT7lcbub89ZgjfVazDuVx30eTMcvreEs4tjbTYlVLO5phyd7tSZpgur5iZPXMKnwbFsaFSC8okRdJ97xp6dNONvpRSWYNjyv3yyP3ChYu8P2UU35erRxzZeTB6GQPDH6VE41YWJ1RKqczjmHJ3u4TmOX6nw/oYtoY0o2LiXjocOkTXzj2tjqaUUpnOL+UuIs2B4YAbGGuMGeqP57nMe/4ca/ZMY949DTEIbaIWMbTVc+Ruks+fT6uUUgErw8tdRNzAp0AT4DCwVkRmGWO2Z/RzAXz11QhGFytBRKnG3Bq/k+eiY3i6nW7Lq5TK2vwxcq8F7DHG7AMQkSlAayDDy/21icP4tnQDgoinXeQC3n/8ZbLnyJHRT6OUUrbjj3IvBURecfswUPvvDxKRrkBXgLJly6btieKSqRq3k9eTs9FIN/pSSqk/+KPcUzv3//9OGzXGjAZGA4SFhaVpO8dXu/RBd1pXSqn/54/tBw4DZa64XRqI8sPzKKWUugp/lPtaIFREyotIMNAWmOWH51FKKXUVGT4tY4xJFJGXgPmkLIUcZ4zZltHPo5RS6ur8ss7dGDMHmOOPn62UUur6HLPlr1JKqT9puSullANpuSullANpuSullAOJMWk6fyhjQ4hEAwf/drgwcNKCOBnJ7q9B81vP7q9B8/tXOWNMkdTuCIhyT42IrDPGhFmdIz3s/ho0v/Xs/ho0v3V0WkYppRxIy10ppRwokMt9tNUBMoDdX4Pmt57dX4Pmt0jAzrkrpZRKu0AeuSullEojLXellHKggCt3ERkgIkdEZJPvo+UV970pIntEZJeINLMy5/WIyOsiYkSksO+2iMgnvvybRaSG1RmvRkQG+jJuEpEFIlLSd9wWr0FEPhSRnb6MM0Qk/xX3BfzvkIg8LiLbRCRZRML+dl/A579MRJr7cu4RkT5W57keERknIidEZOsVxwqKyM8iEuH7XMDKjDfFGBNQH8AA4PVUjt8G/A5kA8oDewG31Xmv8hrKkLLl8UGgsO9YS2AuKVeqqgOstjrnNfLnveLrV4CRdnoNQFPA4/v6A+ADO/0OAbcCVYClQNgVx22R35fV7ctXAQj25b7N6lzXyXwvUAPYesWxYUAf39d9Lv8u2eEj4Ebu19AamGKMiTPG7Af2kHIx7kD0MdCLv15esDUw0aRYBeQXkRKWpLsOY8y5K27m4s/XYYvXYIxZYIxJ9N1cRcrVwMAmv0PGmB3GmF2p3GWL/D61gD3GmH3GmHhgCin5A5Yx5hfg9N8OtwYm+L6eADycqaHSIVDL/SXfn9TjrvgzKLULb5fK/GjXJiIPAUeMMb//7S5b5L9MRAaJSCTwNPCO77CtXoPPs6T8tQH2zH8lO+W3U9ZrKWaMOQrg+1zU4jw3zC8X67geEVkIFE/lrr7A58BAUkaLA4F/kfIP9IYuvJ0ZrpP/LVKmBf7v21I5Ztk61Gu9BmPMTGNMX6CviLwJvAT0J4Bew/Xy+x7TF0gEvr78bak8PmDzp/ZtqRwL1LXMdsrqSJaUuzGm8Y08TkTGAD/5bgbMhbevll9EqpEyF/q7iEBKxg0iUosAyg83/v8A+AaYTUq5B8xruF5+EekAPAg0Mr4JU2yU/yoCJv8NsFPWazkuIiWMMUd9U5AnrA50owJuWuZvc7iPAJffuZ4FtBWRbCJSHggF1mR2vmsxxmwxxhQ1xoQYY0JI+QWvYYw5Rkr+9r4VJ3WAmMt/7gUaEQm94uZDwE7f17Z4DSLSHOgNPGSM8V5xV8D/Dl2HnfKvBUJFpLyIBANtSclvN7OADr6vOwBX+6sq4Fgycr+OYSJSnZQ/4Q4A3QCMMdtE5FtgOyl/ar9ojEmyLOXNm0PKapM9gBfoZG2caxoqIlWAZFJW/HT3HbfLa/gvKStKfvb9BbXKGNPdLr9DIvIIMAIoAswWkU3GmGZ2yQ9gjEkUkZdIWTXmBsYZY7ZZHOuaRGQy0AAoLCKHSflrdSjwrYh0Bg4Bj1uX8Obo9gNKKeVAATcto5RSKv203JVSyoG03JVSyoG03JVSyoG03JVSyoG03JVSyoG03JVSyoH+B60VXY98zUY3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(T):\n",
    "    plt.plot(A_grid, A_policy[i, :])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAD4CAYAAADGmmByAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAaGElEQVR4nO3dbYxc1Z3n8e+/qv0ABDBg82RD7ARvNiSMEGkRtFqNsoEBw0aY7MDK2UhYCZKTEWhfRRMYpCQagjYkO0LKTJIVWdCYKLMOYjeCTWCJycNGswoJZkIIJiF0IECDQwzm2WC7qv77ok7jqu7q7rouN+2yvx+pxK1zz7333HZzf33uObcqMhNJkvZVbb4bIEkabgaJJGkgBokkaSAGiSRpIAaJJGkgI/PdgLfb0qVLc+XKlfPdDEkaKg888MDzmbms17pDLkhWrlzJli1b5rsZkjRUIuLJ6dZ5a0uSNBCDRJI0EINEkjQQg0SSNBCDRJI0kIMiSCJiTUQ8GhFjEXH1fLdHkg4lQx8kEVEHvgZcCJwOfCwiTp/fVknSoeNgeI7kbGAsMx8HiIhNwFrgkXltVYc39+ziuef+yHPPPs32F57jpVdfZuebr7Or0WB3s0GDFk2SBkmLpEXQDMgIEtqvifcBLdrL0FEnIBModQCSAHJvXYCYaNXeem+t61HeacoXDnRtH9PX61G/e7/TrOhH9LftIF+WMCdftDDAKb8dBvo30QFr1Quv8LlPXbvf93swBMly4OmO9+PABzsrRMQGYAPAqaeeul8P/vttT/Kd/3cPT7d28fyiRby06B28Vj+MnbXD2BmH8yaL2RMLS+13wOJ3wOL92gRJB4nI1pzu/y9qP5uT/R4MQdLrT6euPyIz8ybgJoDR0dGB/8B8Y/currz7Vh464iTG6yvguLMBOCJf49jWDo5svs5xe17isMZuFjcaLGw2WdhosqDRYEErqTdbjDSTeraoJ9RaSZ1oLyfUa0Etg3oEtQhqBLVajXqt/ta6qI+0l0t5LdrLtXqdGnWiVmOkDlGrEdSp1YKIoFYfoV5v39Gs1erUo97+Idbrb5XXazVqtfKrUQtGRmrlB12jPrKg64ce9b13R+v1kY7lvf8stXIMgPrIQnoZKfudrD6yd9vp/kaebp9T6tU621qfoeZU0bFtrVZt2+nUat13lhct9i8MzbWz5mSvB0OQjAOndLxfATw7lwf8L//rm9x1wr/lXzV+x9rtP2blk8+z5OVXOXb1uxj98GWsXH4y9Zq3BiQdGg6GILkfWB0Rq4BngHXAf5rLAz5y+CIOz9e59P/8d3538n/k45/8FKced/hcHlKSDlhDHySZ2YiIq4B7gDpwS2Zunctjjh2+gnfv+QM7zriC//rxS1m8YP/c6pCkYTT0QQKQmXcBd70dx3r8hWf4Y/0kznnht1zyb9YaIpIOeUP/HMnb7X8+8FMAVj39IixeMs+tkaT5Z5BU9MDuN1icOzlm7AlGKs78kaSDkUFS0e8Wn8hpjSfYfdhjzsySJAySSp55dQfP1k/m3a88R3PJbkYMEkkySKr4379pPxX6rvHXyJEjGKn745Mkr4QV/PMLL7Ao3+DEra9S37PEHokkcZBM/327/GbBUt7deJyX33iKZuNYx0gkCXskffvTG6/xTP1kTnvtj7xy5FO8mEfbI5EkDJK+3f37BwB497bX2X3s8+zIo+yRSBIGSd+OGv89f/nGbSz/TbJwSYMdHMlIzR+fJDlG0qczGu/lmP+7gG3bfsmy0xrtHkndHokk+Sd1n3YdtYBfbnuA3a1tLN/T4IU80jESScIg6duCxctZeMRHeOEd46xoNHiRIx0jkSQMkr6tPLnFn//zZ3jq2B0sbS2iwQj1Pr8vXJIOZgZJn3aNPcZIcxfPHV/j8PoSagE1eySSZJD0a8EJJ/DrD51Kc+kIO0eWOGNLkgqvhn1a/N73cttFR7J0AewcWeL4iCQVBkmfMpNnXnuGFbt38Xrdp9olaYJB0qdXdr/Ca3teY/kbr/NafQkjPkMiSYBB0rfx18YBWLFnF6/Wj6buGIkkAQZJ38ZfbQfJ8kaDV2t+hLwkTTBI+nTmsjO54fQNnLqnwSs1P7BRkib4WVt9OuGIE7joqNMgk1dqRztGIkmFPZIqXn8egJdrR9sjkaTCIKliZztIXuIox0gkqTBIqnj9eRhZzM5c5KwtSSq8GlaxcwccvpRmpj0SSSoMkip2Pg9HHEejlY6RSFJhkFTx+vPtHknLHokkTTBIqtj5PBxuj0SSOhkkVezcAUeUHonPkUgSMGCQRMRlEbE1IloRMTpp3TURMRYRj0bEBR3la0rZWERc3VG+KiJ+HhGPRcR3ImJhKV9U3o+V9StnO8ac2PMm7H6to0diBksSDN4jeRj4D8BPOwsj4nRgHfA+YA3w9YioR0Qd+BpwIXA68LFSF+AG4MbMXA28CFxRyq8AXszM04AbS71pjzHg+UyvPEPS7pG0WOCtLUkCBgySzPxNZj7aY9VaYFNm7srMJ4Ax4OzyGsvMxzNzN7AJWBsRAXwYuL1svxG4pGNfG8vy7cC5pf50x5gb5al2Dl9Ko+kYiSRNmKv7M8uBpzvej5ey6cqPA17KzMak8q59lfUvl/rT7WuKiNgQEVsiYsv27dv37YwmeiSHH+cYiSR1mDVIIuLeiHi4x2vtTJv1KMt9KN+XfU0tzLwpM0czc3TZsmW9qszujZfa/y2D7Y6RSFLbrJ/+m5nn7cN+x4FTOt6vAJ4ty73KnweWRMRI6XV01p/Y13hEjABHAztmOcb+d8al8K//PdQX0miN+xyJJBVz9Wf1ncC6MuNqFbAa+AVwP7C6zNBaSHuw/M7MTODHwKVl+/XAHR37Wl+WLwV+VOpPd4y5s+AwqNVLj8QgkSQY8PtIIuKjwN8Dy4DvR8SDmXlBZm6NiNuAR4AGcGVmNss2VwH3AHXglszcWnb3WWBTRHwR+CVwcym/GfhWRIzR7omsA5jpGHOt0WrZI5GkYqAgyczvAt+dZt31wPU9yu8C7upR/jg9Zl1l5pvAZVWOMdectSVJezlivA8aftaWJL3FINkHztqSpL28Gu6DRqvlcySSVBgk+8BZW5K0l0GyDxwjkaS9DJKKWq0kE3skklQYJBU1Wu1PYbFHIkltBklFzYkgqfujkyQwSCprtFqAPRJJmmCQVDTRI3GMRJLaDJKKHCORpG4GSUV7eyT+6CQJDJLK7JFIUjeDpKJm0zESSepkkFS0Z2LWlp+1JUmAQVKZs7YkqZtBUlGj6RiJJHUySCpy1pYkdfNqWJFPtktSN4OkIsdIJKmbQVKRz5FIUjeDpCI//VeSunk1rKjhrS1J6mKQVNR0sF2SuhgkFTX8iBRJ6mKQVLR3jMQgkSQwSCpz1pYkdTNIKvLJdknq5tWwInskktTNIKmo0WzP2nKwXZLaBgqSiPhKRPw2Ih6KiO9GxJKOdddExFhEPBoRF3SUryllYxFxdUf5qoj4eUQ8FhHfiYiFpXxReT9W1q+c7RhzyR6JJHUbtEeyGXh/Zv4Z8DvgGoCIOB1YB7wPWAN8PSLqEVEHvgZcCJwOfKzUBbgBuDEzVwMvAleU8iuAFzPzNODGUm/aYwx4PrPys7YkqdtAQZKZP8jMRnl7H7CiLK8FNmXmrsx8AhgDzi6vscx8PDN3A5uAtRERwIeB28v2G4FLOva1sSzfDpxb6k93jDm1t0fiXUFJgv07RvJJ4O6yvBx4umPdeCmbrvw44KWOUJoo79pXWf9yqT/dvqaIiA0RsSUitmzfvn2fTm7CxJPtdZ8jkSQARmarEBH3Aif2WHVtZt5R6lwLNIBvT2zWo37SO7hyhvoz7WumbboLM28CbgIYHR3tWadfjpFIUrdZgyQzz5tpfUSsBz4CnJuZExfpceCUjmorgGfLcq/y54ElETFSeh2d9Sf2NR4RI8DRwI5ZjjFnmn7VriR1GXTW1hrgs8DFmbmzY9WdwLoy42oVsBr4BXA/sLrM0FpIe7D8zhJAPwYuLduvB+7o2Nf6snwp8KNSf7pjzCk//VeSus3aI5nFPwCLgM3t8W/uy8xPZ+bWiLgNeIT2La8rM7MJEBFXAfcAdeCWzNxa9vVZYFNEfBH4JXBzKb8Z+FZEjNHuiawDmOkYc6nZSuq1oJyvJB3yBgqSMiV3unXXA9f3KL8LuKtH+eP0mHWVmW8Cl1U5xlxqlCCRJLU5h7WiZqvl+IgkdTBIKrJHIkndDJKKmq20RyJJHQySito9En9skjTBK2JFjaZjJJLUySCpyDESSepmkFTUbKXf1y5JHQySiuyRSFI3g6SiZtNZW5LUySCpyFlbktTNK2JFzVaLBY6RSNJbDJKKHCORpG4GSUU+2S5J3QySiuyRSFI3g6Sido/EH5skTfCKWJE9EknqZpBU5PeRSFI3g6SiRtMeiSR1Mkgq8rO2JKmbQVKRT7ZLUjeviBU1HCORpC4GSUVNx0gkqYtBUlHDJ9slqYtBUlHT50gkqYtBUpE9EknqZpBU1J7+649NkiZ4RazIWVuS1M0gqcgxEknqZpBU5BiJJHUzSCpotZJMfLJdkjp4Rayg0UoAP2tLkjoMFCQRcV1EPBQRD0bEDyLi5FIeEfHViBgr68/q2GZ9RDxWXus7yj8QEb8u23w1IqKUHxsRm0v9zRFxzGzHmCvNEiSOkUjSXoP2SL6SmX+WmWcC3wM+V8ovBFaX1wbgG9AOBeDzwAeBs4HPTwRDqbOhY7s1pfxq4IeZuRr4YXk/7THmUqPVAnCMRJI6DBQkmflKx9sjgCzLa4Fbs+0+YElEnARcAGzOzB2Z+SKwGVhT1h2VmT/LzARuBS7p2NfGsrxxUnmvY8wZeySSNNXIoDuIiOuBy4GXgX9XipcDT3dUGy9lM5WP9ygHOCEztwFk5raIOH6WY2zr0cYNtHstnHrqqdVOsMOeZhkjMUgk6S2z9kgi4t6IeLjHay1AZl6bmacA3waumtisx65yH8pnbFq/22TmTZk5mpmjy5Ytm2W309vbI3GOgiRNmLVHkpnn9bmvfwK+T3sMZBw4pWPdCuDZUv6hSeU/KeUretQHeC4iTiq9kZOAP5Xy6Y4xZxwjkaSpBp21tbrj7cXAb8vyncDlZWbVOcDL5fbUPcD5EXFMGWQ/H7inrHs1Is4ps7UuB+7o2NfE7K71k8p7HWPOOEYiSVMNOkbypYh4D9ACngQ+XcrvAi4CxoCdwCcAMnNHRFwH3F/q/W1m7ijLfwX8I3AYcHd5AXwJuC0irgCeAi6b6RhzyedIJGmqgYIkM/9ymvIErpxm3S3ALT3KtwDv71H+AnBulWPMlYkeyYhjJJL0Fq+IFTSa3tqSpMkMkgr29kgMEkmaYJBUMDFrq+4YiSS9xSCpwB6JJE1lkFTQcPqvJE1hkFTgrC1JmsorYgX2SCRpKoOkgqYfkSJJUxgkFfgciSRNZZBU4EekSNJUBkkFDaf/StIUBkkFE2Mkfh+JJO3lFbGCht+QKElTGCQVNB0jkaQpDJIKfI5EkqYySCrwyXZJmsorYgX2SCRpKoOkAp9sl6SpDJIK7JFI0lQGSQVNp/9K0hQGSQX2SCRpKoOkgmYrqdeCCINEkiYYJBU0SpBIkvYySCpoNFuOj0jSJAZJBfZIJGkqg6SCZivtkUjSJAZJBe0eiT8ySerkVbGCZqvFAj/5V5K6GCQVOEYiSVMZJBU4RiJJU+2XIImIz0RERsTS8j4i4qsRMRYRD0XEWR1110fEY+W1vqP8AxHx67LNV6M89RcRx0bE5lJ/c0QcM9sx5oo9EkmaauAgiYhTgL8AnuoovhBYXV4bgG+UuscCnwc+CJwNfH4iGEqdDR3brSnlVwM/zMzVwA/L+2mPMZeazfS7SCRpkv1xVbwR+GsgO8rWArdm233Akog4CbgA2JyZOzLzRWAzsKasOyozf5aZCdwKXNKxr41leeOk8l7HmDP2SCRpqoGCJCIuBp7JzF9NWrUceLrj/Xgpm6l8vEc5wAmZuQ2g/Pf4WY7Rq50bImJLRGzZvn17n2c3VbPV8vvaJWmSkdkqRMS9wIk9Vl0L/A1wfq/NepTlPpTP2LR+t8nMm4CbAEZHR2fb77TskUjSVLMGSWae16s8Is4AVgG/KuPiK4B/iYizafcOTumovgJ4tpR/aFL5T0r5ih71AZ6LiJMyc1u5dfWnUj7dMeaMs7Ykaap9vrWVmb/OzOMzc2VmrqR9YT8rM/8I3AlcXmZWnQO8XG5L3QOcHxHHlEH284F7yrpXI+KcMlvrcuCOcqg7gYnZXesnlfc6xpyxRyJJU83aI9lHdwEXAWPATuATAJm5IyKuA+4v9f42M3eU5b8C/hE4DLi7vAC+BNwWEVfQnhl22UzHmEvNVrJwQX2uDyNJQ2W/BUnplUwsJ3DlNPVuAW7pUb4FeH+P8heAc3uUT3uMudJotqgvmqvslaTh5EMRFTQcI5GkKQySCpqOkUjSFAZJBY1W+hyJJE1ikFTQnv7rj0ySOnlVrKDR8jvbJWkyg6SCZtMxEkmazCCpwDESSZrKIKnAWVuSNJVBUkHDwXZJmsKrYgX2SCRpKoOkAmdtSdJUBkkF9kgkaSqDpAI/a0uSpjJI+tRqJZlQd7Bdkrp4VezTnlYLwOdIJGkSg6RPzVb7q94dI5GkbgZJnxolSBwjkaRuBkmfmk2DRJJ6MUj6NNEjqdf9kUlSJ6+KfWp6a0uSejJI+tQos7YcbJekbgZJn+yRSFJvBkmfGk7/laSeDJI+7e2R+COTpE5eFfvUaNojkaReDJI+OUYiSb0ZJH16a9aWn7UlSV0Mkj7ZI5Gk3gySPjlrS5J6M0j61Gg6a0uSevGq2CefbJek3gYKkoj4QkQ8ExEPltdFHeuuiYixiHg0Ii7oKF9TysYi4uqO8lUR8fOIeCwivhMRC0v5ovJ+rKxfOdsx5sLEGMkCB9slqcv+6JHcmJlnltddABFxOrAOeB+wBvh6RNQjog58DbgQOB34WKkLcEPZ12rgReCKUn4F8GJmngbcWOpNe4z9cD49OUYiSb3N1a2ttcCmzNyVmU8AY8DZ5TWWmY9n5m5gE7A2IgL4MHB72X4jcEnHvjaW5duBc0v96Y4xJ3yyXZJ62x9Xxasi4qGIuCUijilly4GnO+qMl7Lpyo8DXsrMxqTyrn2V9S+X+tPta4qI2BARWyJiy/bt2/fpJE84ahEXnXEiRx02sk/bS9LBatarYkTcC5zYY9W1wDeA64As//074JNAr/s/Se/gyhnqM8O6mbbpLsy8CbgJYHR0tGed2XzgncfygXceuy+bStJBbdYgyczz+tlRRHwT+F55Ow6c0rF6BfBsWe5V/jywJCJGSq+js/7EvsYjYgQ4GtgxyzEkSW+TQWdtndTx9qPAw2X5TmBdmXG1ClgN/AK4H1hdZmgtpD1YfmdmJvBj4NKy/Xrgjo59rS/LlwI/KvWnO4Yk6W006A3/L0fEmbRvKf0B+BRAZm6NiNuAR4AGcGVmNgEi4irgHqAO3JKZW8u+PgtsiogvAr8Ebi7lNwPfiogx2j2RdbMdQ5L09on2H/eHjtHR0dyyZct8N0OShkpEPJCZo73WOZdVkjQQg0SSNBCDRJI0EINEkjSQQ26wPSK2A09OKl5K+1mWYTbs5zDs7YfhPwfbP/8O5HN4Z2Yu67XikAuSXiJiy3SzEYbFsJ/DsLcfhv8cbP/8G9Zz8NaWJGkgBokkaSAGSdtN892A/WDYz2HY2w/Dfw62f/4N5Tk4RiJJGog9EknSQAwSSdJADukgiYgvRMQzEfFgeV3Use6aiBiLiEcj4oL5bOdsIuIzEZERsbS8j4j4amn/QxFx1ny3cToRcV1p44MR8YOIOLmUD8U5RMRXIuK3pY3fjYglHeuG4ncoIi6LiK0R0YqI0UnrhuUc1pQ2jkXE1fPdntmUb5T9U0Q83FF2bERsjojHyn+PmWkfB5TMPGRfwBeAz/QoPx34FbAIWAX8HqjPd3unOYdTaH8s/5PA0lJ2EXA37W+RPAf4+Xy3c4b2H9Wx/J+B/zZM5wCcD4yU5RuAG4bwd+i9wHuAnwCjHeVDcQ60v5Li98C7gIWlzafPd7tmafOfA2cBD3eUfRm4uixfPfG7NAyvQ7pHMoO1wKbM3JWZTwBjwNnz3Kbp3Aj8Nd1fM7wWuDXb7qP97ZMn9dx6nmXmKx1vj2DveQzFOWTmD7L9rZ4A99H+pk4Yot+hzPxNZj7aY9WwnMPZwFhmPp6Zu4FNtNt+wMrMn9L+fqVOa4GNZXkjcMnb2qgBGCRwVbktcUtHV3I58HRHnfFSdkCJiIuBZzLzV5NWDUX7J0TE9RHxNPBx4HOleKjOofgk7V4UDGf7JxuWcxiWds7mhMzcBlD+e/w8t6dvg35D4gEvIu4FTuyx6lrgG8B1tP8Kvg74O9oXg+hRf17mSc/S/r+hfWtlymY9yuZtnvdM55CZd2TmtcC1EXENcBXweQ6gc5it/aXOtbS/qfPbE5v1qH/A/htMt1mPsgPxeYFhaedB66APksw8r596EfFN4Hvl7TjtsYcJK4Bn93PT+jJd+yPiDNr3rX8VEdBu479ExNkcQO2H/v8NgH8Cvk87SA6Yc5it/RGxHvgIcG6WG9wcQO2HSv8GnQ6oc5jBsLRzNs9FxEmZua3cxv3TfDeoX4f0ra1J99w/CkzMoLgTWBcRiyJiFbAa+MXb3b6ZZOavM/P4zFyZmStp/890Vmb+kXb7Ly8zn84BXp7oMh9oImJ1x9uLgd+W5aE4h4hYA3wWuDgzd3asOuB/h/owLOdwP7A6IlZFxEJgHe22D5s7gfVleT0wXU/xgHPQ90hm8eWIOJN2N/gPwKcAMnNrRNwGPEL7dsWVmdmct1ZWdxftWU9jwE7gE/PbnBl9KSLeA7Rozzz7dCkflnP4B9qzmjaXnuF9mfnpYfodioiPAn8PLAO+HxEPZuYFw3IOmdmIiKtoz16sA7dk5tZ5btaMIuJ/AB8ClkbEOO1e+JeA2yLiCuAp4LL5a2E1fkSKJGkgh/StLUnS4AwSSdJADBJJ0kAMEknSQAwSSdJADBJJ0kAMEknSQP4/1F5qegPN5bwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(T):\n",
    "    plt.plot(A_grid, V_star[i, :])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To calculate invariant distribution of $A$, we fill in the transition matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "P = np.zeros((N * T, N * T))\n",
    "# P is a transition matrix from (A, Y) to (A', Y)\n",
    "# Rows are current (A, Y), # columns are probabilities of (A', Y)\n",
    "c = 0  # Current row to start from\n",
    "\n",
    "for y_i, y in enumerate(y_grid):\n",
    "\n",
    "    for a_i, A in enumerate(A_grid):\n",
    "        \n",
    "        # Optimal A' given (a_i, y_i)\n",
    "        A_prime = int(A_policy[y_i, a_i])\n",
    "        \n",
    "        P[c, A_prime::N] = π[y_i]  # Probabilities of transitioning to (A', y_i)\n",
    "        \n",
    "        c += 1  # Update current row"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "λ, v = np.linalg.eig(P.T)              # Find eigenvalues\n",
    "unit_λ = np.argmax(np.isclose(λ, 1))   # Locate unit eigenvalue\n",
    "π_star = np.real(v[:, unit_λ])         # Invariant distribution of (A, y) pairs\n",
    "π_star = π_star / π_star.sum()         # Normalise to sum to 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "π_A = π_star.reshape((N, T))\n",
    "π_A = π_A.sum(1)                       # Sum over y states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2deZwdZZnvv885p5ck3VkgAbJBAoQlAQGNEUdFxgVQR9CZwUF0Lq7IvaJzdRzFq6MI3jvuXkVn3OCOoggMjBgBZREEcQgmbCELkBCydBbS2buTdJ+lnvtHVZ2uc/p0+iTp877dXc/388kn51TVqfet03V+9dTvfd6nRFUxDMMwRi8Z3x0wDMMwGosJvWEYxijHhN4wDGOUY0JvGIYxyjGhNwzDGOWY0BuGYYxyTOhTgogcKyLdIpL11P4fRORD0ev3iMi9Q7jv5SJybvT6ahH5+RDu+3+JyE+Gan8H0e47RWRD9Dc7q8Z6FZETD2G/s6LP5gZYP+D3JyLnikjHwbYZffZ9IvLIAdY37PwwTOidIyJrReRNrttV1fWq2qaqpcPdl4j8u4h8+TD68gtVPW+o2lHVear6h0PtT6K9fkKmqv9HVT90uPs+BL4BXBn9zZ700L436j0/jPoxoU8BA0VvI53RelwRxwHLfXfCGB2Y0Hskvp0VkW+IyE4ReVFE3hKtu0REllRt/wkRWRi9fpuIPCkie6Jb/KsT28W35x8UkfXAA9W37CLyfhFZKSJdIrJGRD6S+Py5ItIhIv8oIltFZLOIvD9adznwHuDTka3wmwGO7c0i8qyI7BaR7wFSfdzRaxGRb0ft7BaRpSJy2kDtRHdEnxGRpcBeEcnVuEtqFZFbomN7QkTOSLRdYXnEdw0iMg74LTAtaq9bRKZVWxkicmFkFe2K7IZTE+vWisinomPYHfWhdYDvJyMinxeRddGx/0xEJohIi4h0A1ngaRF5odbnI94a/e22icjXRSRzoH0P0I/ZIvJQ9F3dB0w+QHvxZ/qdF9HyCVFbnVHbn4/7VGMfdZ0f0XsVkStEZJWEv5Pvi4hE67Ii8s3oO3hRRK6sOs/fF31HXdH69wx2fKMSVbV/Dv8Ba4E3Ra/fBxSADxP+sP87sInwpB8LdAFzEp9dDFwSvT4XOJ3wYv0y4CXgHdG6WYACPwPGAWMSy3LRNm8DTojaej2wD3h5Yt9F4BqgCXhrtH5StP7fgS8f4BgnA3uAv40+/4lofx9KHPcj0evzgceBiVFfTgWmDtRO9P09BcwExtT4Tq+OvtO47U8BLwJN0XoFTkzsr9xGdNwdVe1dDfw8en0SsBd4c7TvTwOrgeZEP/4MTAOOAFYCVwzwHX0g+uzxQBvwn8CNifUV/azxeQUejNo5Fng+8f0OuO8a58GjwLeAFuAcwnPu5wO0Odh58TPg10B71M7zwAdr/M3rPj8Sx3on4TlyLNAJXBCtuwJYAcwAJgH3x8dHeO7vAU6Otp0KzPOtAT7+WUTvn3Wq+mMNvfOfEp6MR6vqPsIfzbsBRGQOcAqwEEBV/6Cqz6hqoKpLgV8SCnaSq1V1r6rur25UVe9S1Rc05CHgXuB1iU0KwDWqWlDVu4Fu4OQ6j+mtwApVvU1VC8D/BbYMsG2BUBhOAURVV6rq5kH2/11V3VDruCIeT7T9LaAVOLvOvh+IvwPuUtX7on1/g/Ai+hdVfdukqjuA3wBnDrCv9wDfUtU1qtoNfBa4RA7Ojvqqqu5Q1fWE3/G7D2bfInIs8Ergn1W1V1Ufjvp8IGqeFxIO8v8d8FlV7VLVtcA3gb+vsY+DOT9ivqKqu6JjfZC+7/VdwHdUtUNVdwJfqfpcAJwmImNUdbOqptIOM6H3T/kEj8QdwigM4Cb6fryXAnfE24jIq0Tkweg2eTdhZFN9271hoEZF5C0iskhEdojILsIfX/Lz21W1mHi/L9GvwZiWbFvDcKpmX1T1AeB7wPeBl0TkRyIyfpD9D3hc1etVNQA6oj4dLtOAdVX73gBMT2yTFKwDfWcV+4pe54CjD6I/ye9hHX3HWO++pwE7VXVv1bYHYqDzYjLQXKPd5HeTbLeu8yPBQN9rxb6q9ruX8OJzBbBZRO4SkVMGaWdUYkI/vLkXmCwiZxIK/k2JdTcRRvczVXUC8AMSPmdEzdKkItIC3E4YkR6tqhOBu2t8fiAGK3m6mdBaiduT5Pt+O1P9rqq+AphHaI/80yDtDNZ+su0M4W39pmjRPkJbLOaYg9jvJsJB0njf8XFtHORzg+6L0JIoElpw9ZL8To+l7xjr3fdmYFI0PpHc9lDYRhjtV7db67s5qPNjEDYT/n1jKvajqveo6psJ75SfBX58iO2MaEzohzFR5HQb8HVCL/a+xOp2YIeq9ojIAsKIv16aCT3ZTqAo4QDwwaSzvUTo/w7EXcA8EfnryC74OJWCWkZEXhndnTQR+t89QJwCOlg7A/GKRNv/E+gFFkXrngIujQbxLqDS7noJOHKggUvgVuBtIvLGqL//GO37vw6hj78EPhENhrYB/we4pSpaHox/EpFJIjIT+AfgloPZt6quA5YAXxKRZhF5LfD2QzgWIuvxVuB/i0i7iBwHfBKolZNf9/lRB7cC/yAi00VkIvCZeIWIHC3h4Pk4wr9TN33nVqowoR/+3AS8CfiPqh/q/wCuEZEu4AuEJ3xdqGoX4Y/rVmAn4UVi4UH06XpgroSZJ3fU2P824GJCv3Q7MAf40wD7Gk8YZe0kvNXfTninMWg7B+DXhLfsOwk94r+OvGAIBfHtwC5CL7u8X1V9llAk10RtVtg9qvoc8F7gOsII9u3A21U1fxB9i7kBuBF4mHCwuAf42EHu49eEA9lPEYrn9Yew70uBVwE7gC8SDqgeKh8jvFivAR4hPHdvqN7oIM+Pwfgx4Z3vUuBJwjvTIqGgZwgvxpsIj+/1hL+b1CGhPWYYhjHyie5Of6Cqxw26cYqwiN4wjBGLiIwRkbdKOJ9iOuFdya9892u4YRG9YRgjFhEZCzxEmJ67n9DC+gdV3eO1Y8MME3rDMIxRjlk3hmEYo5xhVxRq8uTJOmvWLN/dMAzDGFE8/vjj21R1Sq11w07oZ82axZIlSwbf0DAMwygjIgPOajbrxjAMY5RjQm8YhjHKMaE3DMMY5ZjQG4ZhjHJM6A3DMEY5JvSGYRijHBN6wzCMUY4JvZFqVr3UxZ9f3OGl7T+u6mT99n2Db2gYh4kJvZFqvvvAaj73q2e8tP2JW57i+kfWeGnbSBcm9Eaq6S2UyJcCT20H3to20oUJvZFqSoFSLPmp4Fr02LaRLkzojVRTDJRS4EdsSx7bNtKFCb2RakqBUvQktsUg8Na2kS5M6I1UUwwCSoF7nzwIlECxiN5wggm9kWp8RfSl6MluRQ8XGSN9mNAbqcaXRx+3aRG94QITeiPV+Iro4zbNozdcYEJvpJpCyU9EX4zy5y2iN1xgQm+kmlIQUAoUVbeCW47oLY/ecIAJvZFqip68cvPoDZeY0BuppuTJK+/z6C3rxmg8JvRGqomtE+cRvad2jXRiQm+kGn8RfeClXSOdmNAbqcY8eiMNmNAbqaZUjqzdeuWWR2+4xITeSDUW0RtpwITeSDXxYKzrfPZCyc+dhJFOTOiNVOMrsi63axOmDAeY0Bupxlf2i3n0hktM6I3UEteEB/PojdGNCb2RWkqJ+jaWdWOMZkzojdSSjKbdR/RWvdJwR11CLyIXiMhzIrJaRK6qsf6TIrJCRJaKyO9F5LjEustEZFX077Kh7LxhHA7JaNq5Rx9n+1jWjeGAQYVeRLLA94G3AHOBd4vI3KrNngTmq+rLgNuAr0WfPQL4IvAqYAHwRRGZNHTdN4xDJ5nxYh69MZqpJ6JfAKxW1TWqmgduBi5KbqCqD6rqvujtImBG9Pp84D5V3aGqO4H7gAuGpuuGcXgko2nXefTm0RsuqUfopwMbEu87omUD8UHgtwfzWRG5XESWiMiSzs7OOrpkGIdPKfA3GBu3rWpRvdF46hF6qbGs5pkpIu8F5gNfP5jPquqPVHW+qs6fMmVKHV0yjMOn4NGjj2fGhm2bT280lnqEvgOYmXg/A9hUvZGIvAn4HHChqvYezGcNwwcVHr1j68Znxo+RPuoR+sXAHBGZLSLNwCXAwuQGInIW8ENCkd+aWHUPcJ6ITIoGYc+LlhmGdyo8ek8zY320baSP3GAbqGpRRK4kFOgscIOqLheRa4AlqrqQ0KppA/5DRADWq+qFqrpDRK4lvFgAXKOqOxpyJIZxkPjNo/d3N2Gkj0GFHkBV7wburlr2hcTrNx3gszcANxxqBw2jURQ9DsZaRG+4xGbGGqllOMyM9dG2kT5M6I3U4nVmrMe7CSN9mNAbqcVnVO1zVq6RPkzojdSSnA1rWTfGaMaE3kgtxYrMF9eDsebRG+4woTdSy7Dx6C290mgwJvRGajGP3kgLJvRGahk+Hr1l3RiNxYTeSC3DZmasRfRGgzGhN1LLsPHoTeiNBmNCb6SWyqjadT16y7ox3GFCb6QWi+iNtGBCb6SWiqjaaz16G4w1GosJvZFahk1Eb3n0RoMxoTdSS9FjLnuxZB694Q4TeiO1xFF1Sy7jPKIvBUpLLlPRD8NoFCb0RmqJvfGWXMa5T15MCL1F9EajMaE3Uks5om/K+onom7IV/TCMRmFCb6SWONMmjOhde/TJiN6ybozGYkJvpJY4km42j94Y5ZjQG6mlFCjZjNCUyTjPoy8GAS25bLkfhtFITOiN1FKMhD6bES8RfWtTFNFbHr3RYEzojdRSCgKaMkIuK56ybiyiN9xgQm+klkLJX0RfLCktTebRG24woTdSSylQctkMuYw4t09Cjz62bizrxmgsJvRGakl69D4ePNKcszx6ww0m9EZqKQUBuYyQy2ScP86vGChNni4yRvowoTdSi++I3tf4gJE+TOiN1FIKNIroPQzGBkouG7ZtM2ONRmNCb6QWi+iNtGBCb6SWUknJZTLksj7SK4OwbfPoDQeY0BupJbZPshn3Rc1i2yibcV9nx0gfJvSGd3oKJZ5cv9N5u31ZN+Il6yYbe/QeSiAsXrvD7iRSRF1CLyIXiMhzIrJaRK6qsf4cEXlCRIoi8rdV60oi8lT0b+FQddwYPdy5dDN/82//xc69eaftVnj0zidMxRG9e9to7ba9XPyDR3no+a1O2zX8kRtsAxHJAt8H3gx0AItFZKGqrkhsth54H/CpGrvYr6pnDkFfjVHKnv0FAoV9hRKTHLZbjDz6JscevapGg7EZL3V29vQUAOjqKTpt1/DHoEIPLABWq+oaABG5GbgIKAu9qq6N1lmemHHQFKISAPmi29On5CnrJm7LV0Qff9+9jr9vwx/1WDfTgQ2J9x3RsnppFZElIrJIRN5RawMRuTzaZklnZ+dB7NoYDcTCU3Bc86UYBFEuu9sB0bitbEa8ZN3ki2F7rr9vwx/1CL3UWHYwZ+axqjofuBT4vyJyQr+dqf5IVeer6vwpU6YcxK6N0UA+8sfTGdG7z7opX1gtok8N9Qh9BzAz8X4GsKneBlR1U/T/GuAPwFkH0T8jBcQC7z6iVy9ZN/4j+vj7tqybtFCP0C8G5ojIbBFpBi4B6sqeEZFJItISvZ4MvIaEt28YkLRu/M1O9RHRN2X9TNYqj4mYdZMaBhV6VS0CVwL3ACuBW1V1uYhcIyIXAojIK0WkA7gY+KGILI8+fiqwRESeBh4EvlKVrWMYHj16Lc9OdevRh8fZF9G7Pe68p+/b8Ec9WTeo6t3A3VXLvpB4vZjQ0qn+3H8Bpx9mH41Rjq8Is5SYGasKQaBkMrWGpIa+XUhk3Ti+k4nvnEzo04PNjDW8U84CcTw4WCgFYVSdDcXdVVQfC3s2qoXv2qP3ZZUZ/jChN7yT9+jRx1F1/N4F8QUlvJtw79HHg7Gus5wMf5jQG96JI/l8qeS03WI8OzUTR/RuhK9U9uj9VK/0NSZi+MOE3vBOX153yiJ6TzNj855mIhv+MKE3vJP3NBhbjD36jEeP3kOtm4LNjE0dJvSGd3xZCcma8PF7V+3CMJgZa4OxqcGE3vCOz5mxcU34+L2rdgFyWT8eva87KMMfJvSGd/ryuj179I7a951H7+vCavjDhN7wjo8yxaraNzM26zbrpv/MWMu6MRqLCb3hHR9T8mNt9ZF1M1zq0bvOcjL8YUJveMeHlRC35cWjr5gZ66HWTdE8+rRhQm94x4d1U535klzWaPry6DOesm4svTJtmNAb3omFJ+9wULKvJnzfzFhXwldKevRZD4Ox5tGnDhN6wzsFD9ZNtU+eXNZoqmvd+BqMtZmx6cGE3vCOjwizOvMlXOZ2MDZu2+XTrcAmTKURE3rDK6rqRei9RvSRwDZlMmQzQhDVwneFDcamDxN6wyulQNFI4/IO0/2q682Ah4g+kfFTUodCb4OxqcOE3vBK0j7wEtFnk1k3riZM+cv4gcSYiHn0qcGE3vBK0j5w69H3pTiWPXpnJRD8jQ+AefRpxITe8Eoy88PHYKwPjz4WWB91dqCyqJk6tIwMf5jQG15JirvTPPqER9/ky6OvGB9wOCs4cXF1PVnL8IMJveGVCqH3MTM263tmrNu7Cai8oNqAbDowoTe8UvDs0Vc+MzZdHj1YYbO0YEJveCVOqRzTlB0GefQpybopBYxpygLQ6/iB7IYfTOgNr8QDg+Nask7T/XzPjM0IZDxF9PliwLiWUOgt8yYdmNAbXimUhT7ndDDWd62bXBTJu76bCILwgSvjWnKA5dKnBRN6wyux0Ixtznny6KUsuu7y6LUs8K4j+kLQ932DDcamBRN6wyuxddPW4tijLyUyX7Lua93EAp91PFkrtmraIuvG6t2kAxN6wyux8PiK6HMenjBVCoLyxSXn+CKTvIMC8+jTggm94ZU4d35cS5ZCSZ3N1Kw9M9bNhaYQJCP6yDZyJPTJwW8w6yYtmNAbXikPxkYRpisrITk7NSuOI/pSf4/eVURfvrA222BsmjChN7yST2TdgDsroZjw6DMZISN+s25clUAoVH3f5tGnAxN6wyux8IxtjqwERxFmsiY8hILv1KOPBL7JtUdfHhOxPPo0UZfQi8gFIvKciKwWkatqrD9HRJ4QkaKI/G3VustEZFX077Kh6rgxOujz6N2m+yVnpwJOn91aDLQ8COvco6/6vu25selgUKEXkSzwfeAtwFzg3SIyt2qz9cD7gJuqPnsE8EXgVcAC4IsiMunwu22MFvo8erfpfsl6MxAKvss8+ly1R++o7Xz1HZRZN6mgnoh+AbBaVdeoah64GbgouYGqrlXVpUD1WXM+cJ+q7lDVncB9wAVD0G9jlFC2Elx79NURfVac1rrJ9vPoXVk35tGnkXqEfjqwIfG+I1pWD3V9VkQuF5ElIrKks7Ozzl0bo4F80U+Emcy6gSiid1jrpl9E71robWZsqqhH6KXGsnrPyro+q6o/UtX5qjp/ypQpde7aGA0USgFNWaElF1k3jjzjWNSbsn2RtUuPPls9M9Z51o3bwW/DL/UIfQcwM/F+BrCpzv0fzmeNFJAvBjRlM+XsE38Rvdusm1yi3WR/Gk3/wW/LukkD9Qj9YmCOiMwWkWbgEmBhnfu/BzhPRCZFg7DnRcsMAwiFvTmXoTmKrF0JT3xBiSdLuYzoC4kJU1nHjzHMV6VXmkefDgYVelUtAlcSCvRK4FZVXS4i14jIhQAi8koR6QAuBn4oIsujz+4AriW8WCwGromWGQYQCk9TNkNTLjwVXVk3yZrw4MGjz3ry6KtnxprQp4JcPRup6t3A3VXLvpB4vZjQlqn12RuAGw6jj8YoplAKaM5myl65yzz62DaBMKIvOmy7X9aNo7bj77elKXyEogl9OrCZsYZX4sHY2KN3WesmFlmIhN6hR9/kqx599P02RRdX8+jTgQm94ZV8sdqjdxRVJ2rCQ1gu2GU9+uqsG1dt9xaTQi82MzYlmNAbXgkjevfWTbImPISlCPx49G5LIMQRfEsuQ3MuY4OxKcGE3vBK9WBsoeiygmQios+4mxlbquHRu54wVbZuLKJPBSb0hlcKxXAwNrZuvHr0DssvVM+MdfcowYCMhMfbnMvYYGxKMKE3vJIvBTTlxL1HX5V1k3OYR5+8yGQygoi7p1vFE9QAG4xNESb0hlfK6ZU5tzNji6Wg7JOD26ybQinoZxu5fJRgc0LozaNPByb0hlf6SiC4nTBVrLJumrIZLxE9uJ6VG2Y5ATRnLY8+LZjQG14plAKacplyhJv3UBMe3Eb0/QeC3WX8FIpaZd2Y0KcBE3rDK4WS0pzNIBL69C49+mw/j97dQHAuWzkr12VEH9tkYdaNefRpwITe8Epo3fQ9P9XlM2P9RfS1PHo3x91bSgzGWh59ajChN7yS9IybHKb7VXv0vrJuwHFEX+wbjDWPPj2Y0BteyZcq0/3cefRBVUSf8ZJHD26fV1txYc1mrARCSjChN7xSKCUjTIcRfclPRB8EiioV4wNZh3V2CiUbjE0jJvSGV5ITeFzO1CwFfYIHodi68OjLDyXP+sm6SY6JhN+3DcamARN6wxulQAmUhJXgzjOu7dE3vu3qRxjGr509SrCfVWYRfRowoTe8kSywFf+fd1bUrNqjdxPRF6KLia+sm0qrzAZj04IJveGNfFno+/K6XUWYvjz6UslvRF89GGvVK9OBCb3hjVhk+qbkuxOeZE14cFePvs+jr5ys5bIefTKP3jz6dGBCb3gjX23d5NxZCaWaM2MdRPSx0Pvy6Iv9PXpVE/vRjgm94Y14+n2zh3S/6lz2WGwbLXqxF19pG7nL4c+XApqjEgjNWbfPqzX8YUJvAKF3+/k7nmHL7h5nbZYj+pyPCVP9Pfp4eaPbTbYHHmrdJC6s8TJXvNDZzZfvXEFgFxenmNAbAKze2s3PF63n4VWdztqMBSaOLN0WNavKunEU3RZrpFfmsg6zboqV9ejBXWlogPtWvMRPHnmRbd29zto0TOiNiO7eYvh/T9FZm/3TK1179D4jel/VK7XvDirn9vGN0Hd+dfW6O88ME3ojIv4Bdjv8Aears25y7rJuitUzYyPhbXhEXyO90tXMWFWtmDDVUrZu3NkoPgIKw4TeiNjTUwAcC32tCVOuCovVyKOHxkf0xQEmTLmI6GNBj62y8uMbHVo3Ps4zw4TeiIh/eF1OrZtQeCpnxpactF1rZmy8vLHtRhF9tnJ8wMms3FLlHZSPwdiydWMRvVNM6A3Aj3VTnjBVUdQsHVk3TR5y+GuVnADHHn2v+/PMMKE3IuIf3l4f1k350XZui5rVjOgbfKGp5dFnHdW6icdEytVCPXr0Ls8zw4TeiIhvpX1k3STT/YqBNjzHulZN+Jyj9Mpy1k228m7CxYSpfI3vG/xYNxbRu8WE3gD6hN5l2lt1hFkWHkc+eXWtG6DhpYprzYx1VWenPCaSuIMC14Ox5tH7wITeAKC7t1DxvwvKWSA5t1ZCrZrwuYzjiN5L1k0c0WcBT3n0Hs4zw4TeiBguE6ag8RFmrRRHZx79AA8eKToQ2747qL6ZyMnljaZQCugphG1ZHr1b6hJ6EblARJ4TkdUiclWN9S0icku0/jERmRUtnyUi+0XkqejfD4a2+8ZQkfROXVUz7D9hKow0G+0ZHyii9zEz1lVEX11bKP7eXQ3GJgdgzaN3S26wDUQkC3wfeDPQASwWkYWquiKx2QeBnap6oohcAnwV+Lto3QuqeuYQ99sYYmLPtFBSeosBrU3ZhrfZ/8EjUrG8UdSqCZ91ZN3UjOhd5dEX/Q7GJn158+jdUk9EvwBYraprVDUP3AxcVLXNRcBPo9e3AW8UEcEYMXR5iLbK1k2mMsJstJUQ2zNJ6yYWvYbPjC35nxlbbZW58uiT4m4RvVvqEfrpwIbE+45oWc1tVLUI7AaOjNbNFpEnReQhEXldrQZE5HIRWSIiSzo73VVPNPro7ikyua25/NoFhVI4OzWT6XuUYLjcfU141zNjaz3dqtGWWfXM2GbHEX0s7pPbmk3oHVOP0NeKzKvPyIG22Qwcq6pnAZ8EbhKR8f02VP2Rqs5X1flTpkypo0vGUFIsBewvlDhmQivgMqKvLCzmykoYKPMlua7xbVd69ACNDuprPaMX3KVXxpk2x0xotcFYx9Qj9B3AzMT7GcCmgbYRkRwwAdihqr2quh1AVR8HXgBOOtxOG0PL3t6wvszUCWMAd/5pvtj3oGpw79HXjuj9ZN2E6xp73Plqj97xYGx8Xk2dMMbKFDumHqFfDMwRkdki0gxcAiys2mYhcFn0+m+BB1RVRWRKNJiLiBwPzAHWDE3XjaEirig41XFEnyyZCwkrocERZu2oOvLoG53DP4BHn+xXoxgondW1Rz91Qiv5YkCvowJ2Rh1ZN6paFJErgXuALHCDqi4XkWuAJaq6ELgeuFFEVgM7CC8GAOcA14hIESgBV6jqjkYciHHoxMLeZ924mcwSPu0oMSDqKMIcqN4MOIzos+7bLlSlV8aD4K49+vg829tboiXX+Owuow6hB1DVu4G7q5Z9IfG6B7i4xuduB24/zD4aDSb+AZYjeoeDsU25YeLRZ/09M7Yc0Tf4Apcv16MPv+dMRshl3BWS6+4pks0IU9payu+PGNfspO20YzNjjbKwlz16T9aNO48+yrqpGVV7GB/Iunm6VbVHD/EzANxF9G0tOdpbmwDosjIIzjChN8oe/eS2ZnIZcRbR54taITotObcRfXVN+OS6RrddK+vGmUef67vIuHwGwJ6eQiT0oZFgmTfuMKE3ytZNe2sTba05pxOmfFg3BZ8efXRsiaad3U0UqqqFxq9dDcZ29xRpb83R1hIJvWXeOMOE3ihHVm0t4Y/QpUdfMRjrqMhW7ZrwjmbGBkpTVkhOHHcd0SfHB5qz4jCPPrRu2lpN6F1jQm/Q3VskIzC2OUtbS86ZR1/o59HHZXN9zoxtvHWTbNdl2/mS0pzLVFxkmnIZp1k3ba052qOI3urduMOE3qCrJ4y0RIT2VncRffWEKfd59LUyXxo/GJv058O23dxN5ItBxZgIhBdXVx59aN00WUTvARN6g67oBwihfeNuwlRVCYRokLDREeZAFSST6xrFASP6BgtueAdV2bZLj35PFFCMacqSEaooM5cAABg9SURBVBuMdYkJvUF3b6E8QNbW2uR0MLY61S9e3kh8Zr4Ug6DiTsJl29VWGUQevTPrpkB7a3jn6DKgMEzoDfq8UwgjelfeaXWEGQte4z16f1k3NSP6rJusm3yp0iqD2LppvNDHT5eKA4r21ibz6B1iQm/QHd1SA6FH72giS75YGWGKCM0OhKdU41GCzrJuSuotoh/Qoy823qOPny5VvnNscXeeGSb0BrFH3/cD7CkEzqK86gizOZdp/DNja+TRxy+dRPRZP1k3Na2bnBuPPo7ey+eZw/kahgm9QVjyICn0UPl8z0ZRHdFDWAbB2czYqruJ8ElPozfrplDSilmx4K4EQj+hdzhfwzChN6i0bmKv3oV/WojyupO4yAIp1PDo4/eNr17ZfzDWb0TvZjC2u2zdRNldre7maxgm9KknfrpU/ANsdzg9fcB0vwZ7xrVqwsfvG11BsljqPxjb59E3/sEjtfPoXQh96MfHgUS7RfROMaFPOfHTpdpaKyP6Rgt9ECjFQGnOVtYjb3YwU7NWTXhwE9GXAq0ovRC3C43Pox8466bxg7FdPbUGY03oXWFCn3LiypXVHn2jo63y80v7ecbuPPp+EX0246TWTbbao3dUC7+WdeNqwlQs9OMTAcW+fKnhx2yEmNCnnHLlykR6JTS+Jn0s5j6shFp59PF7JxH9ANZNwz36ovazylxNmCp79NUBhUX1TjChTzn9f4ChV9/oiD62C2pHmO5rwofvpVxGuFEUg6CGR58pr2skA0X0LqpXdveEhfPGNIVWXbvVu3GKCX3K6a72Tss/wMZOZql+UHVMswPhiSPnKr0lmxEnDx4ZMOvGh0fv6MEjcYniuHKmq4DCCDGhH2YsfHoTz3TsdtZen0cf/vDGNmURBwWnyo+16yc8Ljz6MMUxWa4XoojeiUdfPTbgd2ZsvhSg2ti29/QUyucYuAsokvz5xR3ct+IlZ+0NJ0zohxHFUsCnb3ua6x5Y5azNvqdLhT+8TEYY19z4HOfyYGw/z9iNR1+d+QLDIKL3MBgbP76x0W13J2ZfQ98dpMt6N1/73bP88x3LnLU3nDChH0as2baXnkLA8k17nLVZbd3Erxvv0Q88GNvwomal/rNT47Yb75MruazHmbE1ZiJD45/qFVs3Ma49+lKgrNi8hy17etje3eukzeGECf0wYvmm0LLZuGs/u/e5uaXt7i0i0dOlYlzUIYkLafUTnlyGfLHU0LZrVZAEVxG9v5mx+VJQswQCNL40dLJCKrhL441Zu30v+/LhebVis7tAarhgQj+MWLax7wRcvtmNT598ulSMi8ks+VL4o+tX1MzBBJ5aZQjAo0fvYGasqpIvBrTUuIMCGp5LnyyzAe4m5sUs29j3e3J5xzxcMKEfRizftJtjjxgLwApHJ2NXT5HxiUEyiEsVN3owdqD0SjcTpvxF9H48+njftbKcgIZfXPcknmIGMK7ZrdCv2LSH5myGY8a3mtAb/lBVlm/aw+vmTObo8S3OhD75dKkYpx59DSuh4YOxNWrCQ+iVNzrFMax10z9/H2honZ1yOmuNLCdo/HN646dLxWQzwrjmrDPrZtmm3Zx8TDunz5jAik3ustqGCyb0w4QNO/bT1VPktOkTmDdtgrOoo9o7BTfWzUB59C7K5taqCQ+jO6IfcEzEgUdf/XSpGFc16eMg6rTp45k3bTxrtu1lXz5d+fsm9MOEZVGUMW/aeOZOHc/qzm56Co0dlIT+3ilEP0BHEX2tB2EMhY2gqrzzX//EN+99rt+6WjXhIcxnH4qsm8279/Pya+/jT6u31Wy7+iIjIg2/yOTLWU61B2Mb6dFXP10qpq3FTanijbv2s2tfgbnTJjB36nhUYeXmroa3O5wwoR8mLN+0m1xGOOnoduZNG08pUJ5/qfEnY1dP/4i+vSVHd75I0EDh6R1owtQQefTLN+3hyfW7WPj0pn7rGu3R3/n0ZnbszXPn0lpt1x4IbnSdnbLQ1xj8hsZ69OXKldV3jq1NTqyb+O74tGnjmTd9ApC+zBsT+mHCso17OPGoNlqbssybFp6MLuybrt5iuaJgTFtrDlXYNwR3FKVA6a2RLhkLS3UefXM2SzHQw77I3P3MZgDWbd/Huu17K9Y1Ouvmrqjth5/f1m/Gaa2sm7jtocq62Z+v8X0XB76DgsZaN9WVK2PaHZUqXr5xNxmBU44Zz7QJrUwc25Q6n96EfhgQeoi7OS2KNmZMGkN7S66cV99Ialo3Q1iH5DO3L+Vt332kX6Q8oEcfDw4ehuipKnc9s5njp4wD4I+rKi2UWjXhYWgi+g079vHUhl0cP2UcG3ftZ822yotMKeg/aSlueyguMjf/eT0vv/Y+Onbuq1h+oDERaOxgbPXTpWJcPU5w+aYwiBrTnEVEmDt1fOoyb0zohwFbu3rZ1p1n3rTxQFiG4NRp4xueeVP9dKmYoapD8vxLXdz+RAert3bz0PNbK9YVDlACAQ5vpubyTXtYt30fl7/ueGZMGsPDz3dWtd0/8wWirJvDFNvfLguj+WsvOg2AP1a1XesJU2Hbh3+R6SmU+Pb9z7O/UOKmx9ZXrMsPKPRhX3obGNFXP10qxtVg7LJNu8t3yRCOgz27pavhlUqHEyb0Nbhx0Tq+9Jvlzk6EOHKPI3qAuVPHs3JzV0MH6KqfLhXTPkR1SL5z/yrGNmWZ3NbCjY+uq1g3YFGzIfCM71y6mWxGOH/eMbxuzhQefWF7hTXRUyhRI6Afkoj+rqWbOX36BF5z4mRmHTm24m6iUAooBAFZqXU3cfgXmZseW89Le3qZPXkctyzeUGGZxd9ny0DfdwMj+uqnS8W0teTo6mnsDPBt3b28tKe3HEQBzJ02nnwx4IXOvQf45NDRUyhx1e1L+fVTG520V4u6hF5ELhCR50RktYhcVWN9i4jcEq1/TERmJdZ9Nlr+nIicP3Rdbwy/eGwd/3zHMv7fn9by+TuWNbyqH8DyjXsQgVOn9p2M86aNZ3+hxIvbDv9kvHXJBi763iOsqhrcLVeurJF1A4c3mWXl5j3c9cxmPvDa2Vz6qmP5w/OdrN/eZyfEGUVDne6nqtz9zGZec+JkJo1r5pw5k+nqLfL0hl0AvNDZzeK1O5g/64h+nw09+kMXvA079vF0x27e9rKpAOFFZs328kXt9sc7UIVXzJpUs+3DyaPfny/xbw+9wKuPP5IvXTiP7Xvz/G7Zlor1cCCP/vDO83wx4PN3PMPlP1vSb4wgFvr26oAiiugb+RuLLZrKiD4eA2u8NVosBVx505PcvHgDn7jlqYq/iUsGFXoRyQLfB94CzAXeLSJzqzb7ILBTVU8Evg18NfrsXOASYB5wAfCv0f6GJb95ehOfv2MZf3nyFP77uSdw8+INfPPe54ds/7v3FWqe1Ms27Wb2keMqIp74ZDyc7IAgUL72u2f59G1LeWbjbt57/WMVYlv9GMGYQ6lDki8GFXdA37l/Fe0tOT702uO5dMGxZET4xZ/DqH7L7h7+/b/WcvLR7TUizIMrsrW9u5dfP7WxfEFctnEP63fs469OD8X2L06YTEbg4Siy/s79q2jJZbn8nOP77St7kGK7prObXz3ZUa5LFA8Av+30WOgnsy9f4vF1O+ktlrjugdWcOXMi5540pWbbBxPRB4FWpN/+4rF1dHb18ok3n8Rro7uJ+C6qUAr45r3PMaYpywlHjavYz1Dk0e/eV+CyG/7Mzxet576VL3H5jUsq7iYOdJ4FCvuHYNC/q6dQ8xji0gdzExH98ZPH0ZLLNNwaDQLl07cv5f6VL/HZt5zCGTMn8vGbn2TRmu0NbbcWucE3YQGwWlXXAIjIzcBFwIrENhcBV0evbwO+J2HxlIuAm1W1F3hRRFZH+3t0aLrfR0+hxG2Pd6CEUZ0qiIQ5yhmB3kLA3t4i3fkiTZkM48fkGN/aRCYj9BZK7Nhb4HsPruKVxx3Bv77nFbQ2Zdi1L8/3HlxNT6HE7CnjDti+aihO+VIoeC25LGOas2REeGrDThat2cH6Hfs4Znwrbzz1KP7y5KPIlwJWvdTNYy/u4LUnTq7Y34lHtdGUFX795MYD3t4Kwv5Cia1dPWzd00u+GHDckWOZPXkcDz3fyZ1LN/PuBTN579nH8Z6fPMalP1nErR95NY++sJ2v3/McIvQ7tljo71+5lR378jXbzYgghBU3l6zdwbKNexjTnOX8eUdz5sxJ/G75Fv7hjXOYMLaJCTRx3tyjuXXxBj72hjl89KYn2JcvcfPlZ/WrCR9HmP/5xEYmtzfXbDcj4Xd9/8qtPLJ6Wzld8p1nTUc1jI7Pm3c0ABPGNnHmzIn8cVUnbzt9Kr9ZuomPnHMCk9ta+u07lxW6eor84rF1/dbFbWdF2La3l7uWbi5Hi+2tOT78uuP53bItnDFjAjOjMhavPuFIchnhj6s6Wb21i4279vMvf316v2OO217d2X3AtoVwPGfJup08uX4nvYWAc06awltPP4YfPPQCrz1xMgtmh3cq7z37OL5810pWbNrD7U90sGTdTr5zyZlMnTCmYr/xhfXhVZ3sPcAkIkHYvb9QPs9amjIcP3kcMyaN5boHVrF+xz6+9a4zKJZCcfv4L5/kq3/zMn748Bquf+RFpk5oLT9dKia+c7zx0XX97MNqCsWA3mJQLrU8tjlLSy7Lyi17WLRmB89u2UNbS45zTz6KN516FG0tOVZv7eaOJzdy7BFjmTCmbxwql81wyjHt/HHVtgG/7xjVSE+i15lIUyC8S+ruLbK/UGJMU5YJY5poa80RBEpvMeDJ9Tu546lNfPLNJ/GR15/Au+bP5OIfPsqHf7qET553Uj/bEuDIcS1ccNoxB+zToVCP0E8HNiTedwCvGmgbVS2KyG7gyGj5oqrPTq9uQEQuBy4HOPbYY+vtewXdvUU+X0et6eZc+ADoWl7smTMn8pP3zWdMVMnxy+84nT37i/zkkRcPqU8xE8Y0sWD2Ebxr/gyWbdzDr57cyC+iwTKRMMvmnWdVfi3NuQyvOG4Sv392K79/dmut3VZun80wpb2F5lyGe5ZvKUeHn7ngFK54/fGICD/7wAIu/fFjnPv1P5AvBZwxYwLXXXoWpxwzvmJfk9taaG/JcfsTHdz+RMeg7Z4+YwLve80stnX1cvczW7h1SQftrTk+8NrZ5e3+/uzj+O2yLVz8g0dZuXkP1737LE48qr3f/mIh+vb9g99JzZg0hsvPOZ43nnIUv122hRsXrSNfDHj9SVOYOLbvIvG6OVO47oFVfPmuFYxrzvGRGtF83HZXb5HP/Wrw8+iMmRP557+ay7xp47n+kRf51n1hfz/7llPK27S3NvHyYyfx+5Vb2bU/zytnTeJ1cybX3N/UCa0sWrOjbDENhAicfHQ7bz9jGi25DL99Zgv3rwwfpvGJN88pb3fxK2byjXuf4x//42lWbt7D+/5iFhed2e+nx8SxzYxrzvKfT2zkP58Y3ENua8lxVHsL+wul8vbjW3P87AOv4tUnHAnAvnyRq3+zggef/T35UsA7z5rOP51/cr8L3HFHhAHGv/z22UHbHYjWpvB38vE3zGHz7v088OxWfpOYO3FUewvvf83sfp87+/gj+eHDa+r6Ww9GSy5TnhdSzUfOOZ6PveFEACaNa+ZnH1jA3/3oUb70mxU1tz9z5sSGCL0M5o+JyMXA+ar6oej93wMLVPVjiW2WR9t0RO9fIIzcrwEeVdWfR8uvB+5W1dsHam/+/Pm6ZMmSgz6QUqBs7+5FRMJIPl4eRfctuQzjWnI0ZTOoKnvzJXbvD62UllyWlqYMbc05MlUZEarKtu58XT5icy5DSy5LNiP0Fkvsz5fIlwKmTRhTsd+eQomnNuyirSXHCVPayheWavLFgF0DRNQAcY+asxkmjm0q/5CKpYCOnftRYPbkymh98dodfP13z3Hpq47lwjOm9TvemL29xfKMxlrtqkKgyhHjmmlNRGo9hRIPP9/JkW0tvOK4Pi9aVXnjtx5iTede3v+aWXzx7fMGPK5d+/I1rZu43VL0t5g2obVCPLbs7uHni9Zx3ryjedmMieXlj6/bwd/8W3gT+fE3nMgnzzu59nGp0tnd2/fFVrUdaBggtOSyTGmvvCN4esMu7ly6iY/+5YkVF5nrfr+Kb0YXgV9++OyyGFZzoL918vtub81VFAcLAuXx9TvZ3t3LBadNrfjcp297mluXdHDWsRO55fJX14wg4cB/62T77a05xiXsxX35Imu37ePo8S0cWXWHdP0jL/LHVZ184k0nccbMiQzEzr35QW0jJTzHm3OZci2kffkSPYUSR49vrTiuIFCe2bibYhBw4lHtFZF8xT5V6eyqry597ArEfQlUQWFMc5axzbnyIH5XT4GuniLZjNDalKW1KcPY5v6x9IH+1rlshiPG9b+TrbOfj6vq/Jrr6hD6VwNXq+r50fvPAqjqvyS2uSfa5lERyQFbgCnAVcltk9sN1N6hCr0x/Hnwua3cu/wlvnThvAFFpxEUSwFnXXsfAvzxM28Y8MffCJ7asIt3fP9P/MUJR3LTh8921i7A2m17+ca9z/G5t53az7IxRh+HK/Q54HngjcBGYDFwqaouT2zzUeB0Vb1CRC4B/lpV3yUi84CbCKP7acDvgTmqOuDoiwm90QhuXbKB8a25flFvowkC5ev3PsffvHx6TavKMIaKAwn9oB595LlfCdwDZIEbVHW5iFwDLFHVhcD1wI3RYOsOwkwbou1uJRy4LQIfPZDIG0ajeNf8mV7azWSEz1xwyuAbGkYDGTSid41F9IZhGAfPgSJ6mxlrGIYxyjGhNwzDGOWY0BuGYYxyTOgNwzBGOSb0hmEYoxwTesMwjFGOCb1hGMYoZ9jl0YtIJ1CrpNxkYFuN5SMF679/RvoxWP/9M5yP4ThV7V8Dm2Eo9AMhIksGmgwwErD++2ekH4P13z8j9RjMujEMwxjlmNAbhmGMckaS0P/IdwcOE+u/f0b6MVj//TMij2HEePSGYRjGoTGSInrDMAzjEDChNwzDGOUMa6EXkatFZKOIPBX9e2ti3WdFZLWIPCci5/vsZz2IyKdEREVkcvReROS70TEsFZGX++5jLUTk2qh/T4nIvSIyLVo+Uvr/dRF5Nurjr0RkYmLdiDiHRORiEVkuIoGIzK9aN1KO4YKoj6tF5Crf/akHEblBRLaKyLLEsiNE5D4RWRX9P+lA+xg2qOqw/QdcDXyqxvK5wNNACzAbeAHI+u7vAY5jJuETutYBk6NlbwV+S/gc87OBx3z3c4C+j0+8/jjwgxHW//OAXPT6q8BXR9o5BJwKnAz8AZifWD4ijoHwyXQvAMcDzVGf5/ruVx39Pgd4ObAssexrwFXR66vi82m4/xvWEf0BuAi4WVV7VfVFYDXhc2mHK98GPk34EPmYi4CfacgiYKKIuH2gaR2o6p7E23H0HcNI6f+9qlqM3i4CZkSvR8w5pKorVfW5GqtGyjEsAFar6hpVzQM3E/Z9WKOqDxM+GjXJRcBPo9c/Bd7htFOHyEgQ+iuj2+4bErdJ04ENiW06omXDDhG5ENioqk9XrRpJx/C/RWQD8B7gC9HiEdP/BB8gvAuBkdn/akbKMYyUftbD0aq6GSD6/yjP/amLQR8O3mhE5H7gmBqrPgf8G3AtYRR5LfBNwh+r1NjeW57oIMfwvwjtg34fq7HMyzEcqP+q+mtV/RzwORH5LHAl8EVGUP+jbT5H+ID6X8Qfq7H9sDyH4mOo9bEay4ZjvvRI6eeoxbvQq+qb6tlORH4M3Bm97SD0vWNmAJuGuGt1M9AxiMjphN7p0yICYT+fEJEFDKNjqPdvANwE3EUo9COm/yJyGfBXwBs1MlcZRv2Hg/obJBlWx3AARko/6+ElEZmqqpsjq3Kr7w7Vw7C2bqo833cC8ej3QuASEWkRkdnAHODPrvs3GKr6jKoepaqzVHUW4Qn/clXdQngM/y3KXjkb2B3fEg4nRGRO4u2FwLPR65HS/wuAzwAXquq+xKoRcQ4Nwkg5hsXAHBGZLSLNwCWEfR+JLAQui15fBgx0tzWs8B7RD8LXRORMwtu8tcBHAFR1uYjcCqwgvB3/qKqWvPXy0LibMHNlNbAPeL/f7gzIV0TkZCAgzBq6Ilo+Uvr/PcKslPuiu6pFqnrFSDqHROSdwHXAFOAuEXlKVc8fKcegqkURuZIw8ywL3KCqyz13a1BE5JfAucBkEekgvJP9CnCriHwQWA9c7K+H9WMlEAzDMEY5w9q6MQzDMA4fE3rDMIxRjgm9YRjGKMeE3jAMY5RjQm8YhjHKMaE3DMMY5ZjQG4ZhjHL+Py03RWy6TZy9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(A_grid, π_A)\n",
    "plt.title('Invariant distribution of bold holdings')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}