2.1 - Ramsey IRF.ipynb

{
 "metadata": {
  "name": "2.1 - Ramsey IRF"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#### 2.1 - Ramsey IRF"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import scipy.linalg\n",
      "\n",
      "# Model parameters\n",
      "a = 0.5;  # capital's share of income\n",
      "b = 0.95; # discount rate\n",
      "d = 0.1;  # depreciation rate\n",
      "s = 1.01; # coefficient of relative risk aversion\n",
      "\n",
      "# Steady-state values\n",
      "kss = (((1/b) + d - 1)/a)**(1/(a-1));\n",
      "css = kss**a - d*kss;\n",
      "\n",
      "# Functional specifications\n",
      "def u(c):\n",
      "  return (1 / (1 - s)) * (c**(1 - s) - 1);\n",
      "def du(c):\n",
      "  return c**(- s) - (1 / (1 - s));\n",
      "def f(k):\n",
      "  return k**a;\n",
      "def df(k):\n",
      "  return a*(k**(a-1));\n",
      "\n",
      "# Linearization\n",
      "phi = b*a*(1-a)*(kss**a)/s;\n",
      "B = np.matrix([[1,phi],[0,1]]);\n",
      "A = np.matrix([[1,0],[(-css/kss),1 + a*(kss**(a-1)) - d]]);\n",
      "M = np.linalg.inv(B)*A;\n",
      "\n",
      "# Eigenvalue decomposition\n",
      "vals, vecs = np.linalg.eig(M);\n",
      "S = np.linalg.inv(vecs);\n",
      "Lambda = np.diag(vals);\n",
      "\n",
      "# Generalize Schur decomposition (sort not working)\n",
      "#S2, T, Q, Z = scipy.linalg.qz(B, A, sort='ouc');\n",
      "#ree_schur_coef = -Z[1,1] / Z[1,0];\n",
      "\n",
      "# Is it saddle path stable (determinate), with largest eig first?\n",
      "test = np.linalg.norm(M - np.linalg.inv(S)*Lambda*S);\n",
      "if not (test < 10e-8 and vals[0] > 1 and 1 > vals[1] and vals[0] > vals[1]):\n",
      "  display(\"Did not work.\");\n",
      "reecoef = -S[0,1] / S[0,0];\n",
      "\n",
      "#display(reecoef)\n",
      "#display(ree_schur_coef)\n",
      "\n",
      "# Time path in two ways\n",
      "k0 = 0.6;\n",
      "c0 = reecoef*k0;\n",
      "T = 120;\n",
      "\n",
      "M1 = (-css/kss)*reecoef + 1 + a*(kss**(a-1)) - d;\n",
      "\n",
      "# Iterate to get paths under each method\n",
      "# under method 1, where we take powers of the inverse matrix times the initial vector\n",
      "# to get the time paths of both consumption and capital.\n",
      "# note that this is vulnerable to computational rounding error problems\n",
      "x = [np.matrix([c0, k0]).getT()];\n",
      "# under method 2, where we take the linearized solution for ct (in terms of kt) and\n",
      "# plug it directly into the capital accumulation equation, then we can simply iterate\n",
      "# that equation to get the path of capital.\n",
      "k_alt = [k0];\n",
      "\n",
      "# Set iterations\n",
      "T = 120;\n",
      "\n",
      "# Get capital paths\n",
      "for i in range(T-1):\n",
      "  x.append(\n",
      "    np.dot(M,x[-1])\n",
      "  );\n",
      "  k_alt.append(\n",
      "    k_alt[-1]*M1 # ???\n",
      "  );\n",
      "    \n",
      "# Extract the capital path for the first method\n",
      "k = [t[1,0] for t in x];\n",
      "\n",
      "# Plot the results\n",
      "fig   = plt.figure();\n",
      "axes  = fig.add_subplot(1,1,1);\n",
      "\n",
      "# and a line showing equilibrium capital\n",
      "axes.plot(range(T), [0]*T, color='gray');\n",
      "\n",
      "# plot capital path from method 1\n",
      "line, = axes.plot(range(T), k, color='red'); # note the tuple unpacking here\n",
      "# plot capital path from method 2\n",
      "line, = axes.plot(range(T), k_alt, color='blue', linestyle='dashed', linewidth=2); # note the tuple unpacking here\n",
      "\n",
      "# add a stem plot\n",
      "\n",
      "# prettify axes and lables\n",
      "axes.set_title('Ramsey Capital Path');\n",
      "axes.set_xlabel('Time');\n",
      "axes.set_ylabel('Capital');\n",
      "axes.set_ylim(-0.2, 1);\n",
      "\n",
      "# Show the plot\n",
      "#display(fig);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEVCAYAAAAYZ2nCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcVGX+B/DPGS6BggoioAwICnLzAoaZaYaaoa6SZhfM\nvK1LZllZuutvu6zalsq2VpbZmuulfl5yKws3FY0E/ZkKKqaJpmQQIygiFwUFgZnn98fkJHIdGObM\nmfm8X695yZl5OPN9AOczz3Oec0YSQggQEREZQSV3AUREpDwMDyIiMhrDg4iIjMbwICIiozE8iIjI\naAwPIiIyGsODyMJs2rQJMTExLf7+6OhorF271oQV1W/69Ol4/fXX2/x5yDIxPMgs/P390a5dO7i6\nusLb2xtTpkzBtWvX5C6rVTZv3oyoqCi4urqiW7duGDNmDL7//vtW73fy5MnYvXu3YVulUuGXX35p\n9vdLkgRJkup9bNGiRXBwcICrqyvc3NwwePBgHD58uMl9btiwAffff3+zn4esH8ODzEKSJHzzzTco\nKyvDiRMn8OOPP+LNN9+Uu6wWe+edd/DSSy/htddew+XLl6HRaPDcc89h+/btbfJ8pjqXV5IkTJo0\nCWVlZSgsLMSQIUPwyCOPmGTfZFsYHmR2Xl5eeOihh5CZmWm4b9myZQgMDESHDh0QHh6Or7/+2vDY\nhg0bMHjwYLz88stwc3NDYGAgDh48iPXr18PPzw9eXl749NNPDe137tyJ8PBwdOjQAWq1GsuXLzc8\n9s033yAiIsLwrvvHH38EALz99tt49NFHa9X5wgsvYO7cuXXqv3r1KhYuXIhVq1Zh/PjxcHZ2hp2d\nHf7whz8gISEBAJCeno5BgwbBzc0N3bp1w/PPP4/q6mrDPlQqFT744AP07NkTXbp0wV/+8hdDQNz+\nLn/o0KEAgH79+sHV1RWff/45SktLMXbsWHh6esLd3R3jxo1DXl5es372QgjD89jb22Pq1Km4dOkS\nioqKGvwdnDlzBrNnz8ahQ4fg6uoKd3d3w/6Ki4sxduxYdOjQAffee69RIyRSOEFkBv7+/iI5OVkI\nIYRGoxF9+vQRixcvNjz++eefi4sXLwohhNi6dato3769uHTpkhBCiPXr1wt7e3uxYcMGodPpxGuv\nvSZ8fHzEnDlzRFVVldizZ49wdXUV169fF0II4e3tLQ4cOCCEEKK0tFRkZGQIIYTIyMgQnp6eIj09\nXeh0OvHJJ58If39/UVVVJfLz80X79u1FaWmpEEKI6upq4enpafje2+3atUvY29sLrVbbYH+PHTsm\n0tLShFarFTk5OSI0NFS89957hsclSRLDhw8XJSUlIjc3V/Tq1Uv8+9//NvR3yJAhtdqeP3/esF1U\nVCS2bdsmKioqRFlZmXjsscfE+PHjDY9HR0eLtWvX1lvXwoULxVNPPSWEEKKyslLMnz9fdO/evcnf\nwYYNG2rVJIQQ06ZNE507dxZHjhwRNTU1YvLkySIuLq7BnwlZF4YHmUX37t2Fi4uLcHV1FZIkifHj\nxzf64hsRESESExOFEPoX06CgIMNjJ0+eFJIkicuXLxvu69y5szhx4oQQQgg/Pz+xevVqcfXq1Vr7\nfOaZZ8Trr79e677g4GCxf/9+IYQQo0aNEmvWrBFCCPHf//5XhIeH11vbxo0bhbe3d3O7LoQQ4t13\n3xUTJkwwbEuSJHbv3m3YXrVqlRgxYoShv42Fx52OHz8u3NzcDNtNhYejo6Po1KmT8PT0FCNGjKg3\nIIWo+zu4MzymT58u4uPjDds7d+4UISEhDdZJ1oXTVmQWkiQhMTER165dQ2pqKvbu3YujR48aHv/0\n008RGRkJNzc3uLm54dSpUygqKjI87uXlZfja2dkZANClS5da95WXlwMAvvzyS+zcuRP+/v6Ijo42\nHBD+9ddfsXz5csNzuLm54cKFC8jPzwcATJs2DRs3bgQAbNy4EVOmTKm3L507d8aVK1eg0+ka7O+5\nc+cwduxYdO3aFR07dsSrr75aqz8A4Ovra/jaz8/PUEdTbty4gVmzZsHf3x8dO3bEAw88gKtXrzb7\nuMgTTzyBkpISFBQUIDk5GZGRkQCa/h3U587fy63fAVk/hgeZ3dChQ/H8889jwYIFAPQv6k8//TQ+\n/PBDFBcXo6SkBL17927xQeKoqCh8/fXXKCwsxPjx4/H4448D0L9Av/rqqygpKTHcysvL8cQTTwAA\nHn74YZw8eRKnTp3Cjh07MHny5Hr3P2jQINx111346quvGqxh9uzZCAsLw88//4yrV6/irbfeqhM2\nubm5tb728fFpVv+WL1+Oc+fOIT09HVevXsW+fftqHctojCRJ9bZr6nfAVVV0J4YHyWLu3LlIT09H\nWloarl+/DkmS4OHhAZ1Oh/Xr1+PUqVMt2m91dTU2bdqEq1evws7ODq6urrCzswMAxMfH41//+hfS\n09MhhMD169exY8cOw7tlZ2dnTJw4EU8++SQGDhwItVpd73N07NgRb7zxBp577jkkJibixo0bqK6u\nxq5duwyBWF5eDldXV7Rr1w4//fQTPvroozr7+ec//4nS0lJoNBq8//77hhC7k5eXF86fP2/YLi8v\nh7OzMzp27Iji4mIsXry4zvc0FCQN3d/U78DLywsXLlyoddC/peFO1oHhQbLw8PDAtGnTkJCQgLCw\nMMybNw+DBg2Ct7c3Tp06hSFDhhja1nc+QWPvhDdu3IiAgAB07NgRH3/8MTZt2gQAuPvuu7FmzRrM\nmTMH7u7uCAoKqrVKC9BPXZ06darBKatbXn75Zbzzzjt488034enpCT8/P6xatQoTJkwAoA+GzZs3\no0OHDnj66acRFxdXp+aHH34Yd999NyIjIzF27FjMnDmz3v4uWrQI06ZNg5ubG7744gvMnTsXFRUV\n8PDwwH333YfRo0c3++fT0LkZTf0ORowYgfDwcHh7e8PT07PBfXGEYjskwbcPRAYajQYhISEoKCiA\ni4tLmz2PSqXCzz//jB49erTZcxC1JY48iH6j0+mwfPlyTJo0qU2Dg8ga2MtdAJEluH79Ory8vBAQ\nEICkpKQ2fz5O75DScdqKiIiMpuiRB9+9ERG1TGvHDYo/5nFrfbs13hYuXCh7Dewf+6eY/p07BxEY\naJ19M/HNFBQfHkREAIDKSsDJSe4qbAbDg4isA8PDrBgeFiw6OlruEtoU+6dsFtc/E4aHxfXNAil6\ntVVD1+khIhv07bfAP/6h/5caZYrXTo48iMg6cNrKrBgeRGQdGB5mxfAgIuvA8DArhgcRWQeGh1kx\nPIjIOlRWAr99yiS1PYYHEVkHjjzMiuFBRNaB4WFWDA8isg4MD7NieBCRdWB4mBXDg4isQ0UFw8OM\nGB5EZB3KywFXV7mrsBkMDyKyDmVlAD973mwYHkRkHcrLGR5mxPAgIuvA8DArhgcRWQce8zArhgcR\nWQce8zArhgcRWQdOW5kVw4OIrAPDw6wYHkSkfNXVgFYL3HWX3JXYDIYHESnfrVGHJMldic1geBCR\n8nHKyuwYHkSkfAwPs2N4EJHy8RwPs2N4EJHy8RwPs2N4EJHycdrK7BgeRKR8DA+zY3gQkfIxPMyO\n4UFEyldWxgPmZiZbeCQlJSEkJARBQUFISEio8/iVK1cwatQoREREoHfv3tiwYYP5iyQiZeDIw+xk\nCQ+tVos5c+YgKSkJp0+fxpYtW3DmzJlabVauXInIyEj88MMPSE1Nxbx581BTUyNHuURk6RgeZidL\neKSnpyMwMBD+/v5wcHBAXFwcEhMTa7Xp2rUrrl27BgC4du0aOnfuDHt7eznKJSJLx/AwO1lejfPy\n8uDr62vYVqvVSEtLq9UmPj4ew4cPR7du3VBWVob//Oc/9e5r0aJFhq+jo6MRHR3dFiUTkSXjMY9G\npaamIjU11aT7lCU8pGZcvGzJkiWIiIhAamoqzp8/j5EjR+LEiRNwveMP5PbwICIbxZFHo+58Y714\n8eJW71OWaSsfHx9oNBrDtkajgVqtrtXm4MGDeOyxxwAAPXv2REBAAM6ePWvWOolIIRgeZidLeERF\nRSErKws5OTmoqqrC1q1bERsbW6tNSEgIkpOTAQAFBQU4e/YsevToIUe5RGTpGB5mJ8u0lb29PVau\nXImYmBhotVrMnDkToaGhWL16NQBg1qxZeOWVVzBjxgz069cPOp0O//jHP+Du7i5HuURk6XhhRLOT\nhBBC7iJaSpIkKLh8IjKV7t2B/fv1/1KTTPHayTPMiUj5OG1ldgwPIlI+hofZMTyISNmqqgAhgLvu\nkrsSm8LwICJl46hDFgwPIlI2hocsGB5EpGwMD1kwPIhI2XiOhywYHkSkbGVlHHnIgOFBRMrGaStZ\nMDyISNkYHrJgeBCRsvGYhywYHkSkbDzmIQuGBxEpG6etZMHwICJlY3jIguFBRMrG8JAFw4OIlK2s\njAfMZcDwICJl48hDFgwPIlI2hocsGB5EpGwMD1kwPIhI2XjMQxYMDyJSNo48ZMHwICJlY3jIguFB\nRMrG8JAFw4OIlOvmTUCSAEdHuSuxOQwPIlIujjpkw/AgIuVieMiG4UFEysXwkA3Dg4iUi+d4yEa2\n8EhKSkJISAiCgoKQkJBQb5vU1FRERkaid+/eiI6ONm+BRGT5OPKQjb0cT6rVajFnzhwkJyfDx8cH\nAwYMQGxsLEJDQw1tSktL8dxzz2H37t1Qq9W4cuWKHKUSkSVjeMhGlpFHeno6AgMD4e/vDwcHB8TF\nxSExMbFWm82bN2PixIlQq9UAAA8PDzlKJSJLxo+glY0sI4+8vDz4+voattVqNdLS0mq1ycrKQnV1\nNYYNG4aysjK8+OKLmDJlSp19/W3ESOgGDYa9PRAdHc3pLSJbUlICuLnJXYXFS01NRWpqqkn3KUt4\nSJLUZJvq6mpkZGTgu+++w40bNzBo0CDce++9CAoKqtXuzb27sWgo8LeFPPZPZHOKi4HOneWuwuLd\n+cZ68eLFrd6nLOHh4+MDjUZj2NZoNIbpqVt8fX3h4eEBZ2dnODs7Y+jQoThx4kSd8BBQofBCBQBn\nc5RORJakqAi44zWBzEOWt+tRUVHIyspCTk4OqqqqsHXrVsTGxtZq8/DDD+PAgQPQarW4ceMG0tLS\nEBYWVu/+CvOqzVE2EVmaoiKOPGQiy8jD3t4eK1euRExMDLRaLWbOnInQ0FCsXr0aADBr1iyEhIRg\n1KhR6Nu3L1QqFeLj4xsOjwKdOcsnIktRXAy4u8tdhU2ShBBC7iJaSn/sRKBP96s4mdNR7nKIyNwG\nDAA+/BC45x65K1EUSZLQ2pd+qzjKXFWp2PwjotbgyEM2skxbmVLVX16DQwdnAK/KXQoRmRuPechG\n8SMPBy93oLBQ7jKIyNxqavRnmHfklLUcFB8e6NKF4UFki0pKgE6dAJXyX8aUSPk/dYYHkW3ilJWs\nGB5EpEw8WC4rxYeH8OiCsoIbKC+XuxIiMiuOPGSl+PB49cOu6HDxLN5fweW6RDaFIw9ZKT483D0d\nAACF+VUyV0JEZsWRh6wUHx63PuajMI/hQWRTOPKQleLDo0sX/b+Fl7TyFkJE5sWRh6ysJzy44IrI\ntjA8ZGUV4eGgqoGdjpdlJ7IpnLaSleKvbeXvD9x86a+QPLsA+Ivc5RCRuXDkISvFh4ckAfDkiYJE\nNocjD1kpftoKgH7JFcODyLZw5CGrBkcexcXFjX6juyUlPi9RQmRbbt4EqqoAFxe5K7FZDYZH//79\nf/ukvvplZ2e3SUEtwvAgsi23pqwaeY2ittVgeOTk5JixjNapceuCoksCd5Xqr9BMRFauqIjHO2TW\nrAPmJSUlyMrKQmVlpeG+oUOHtllRxvrLB2q8qzmCf6wB/vxnuashojZXXMzjHTJrMjzWrFmD999/\nHxqNBpGRkTh8+DAGDRqEvXv3mqO+ZvHo6ggAKLxYDcBB3mKIqO3xYLnsmlxttWLFCqSnp8Pf3x8p\nKSk4fvw4OlrYxz528dTPexZeuClzJURkFlymK7smw8PJyQnOzs4AgMrKSoSEhODs2bNtXpgxbl2i\n5MqlGnkLISLz4MhDdk1OW6nVapSUlGD8+PEYOXIk3Nzc4O/vb4bSmo/XtyKyMRx5yK7J8Pj6668B\nAIsWLUJ0dDSuXbuGUaNGtXlhxvD0BNwcy+FiVwGAy62IrF5REdCjh9xV2LQmp62mTJli+Do6Ohqx\nsbGYOXNmmxZlrKAgoHj2a0j+4xa5SyEic+BSXdk1GR6nTp2qtV1TU4Njx461WUEtxhMFiWwHl+rK\nrsHwWLJkCVxdXfHjjz/C1dXVcPP09ERsbKw5a2weXt+KyHbwgLnsGgyPV155BWVlZZg/fz7KysoM\nt+LiYixbtqzVT5yUlISQkBAEBQUhISGhwXZHjhyBvb09tm3b1vgOu3QBrlxpdV1EpAA8YC67Bg+Y\n//TTTwgJCcFjjz2GjIyMOo/379+/xU+q1WoxZ84cJCcnw8fHBwMGDEBsbCxCQ0PrtFuwYAFGjRoF\nIUTjO+3aFcjPb3FNRKQQQnDkYQEaDI/ly5djzZo1mDdvXr0XSExJSWnxk6anpyMwMNCw5DcuLg6J\niYl1wuODDz7Ao48+iiNHjjS5zypvP1zIluCUD3Tr1uLSiMjSVVToL4j42/lnJI8Gw2PNmjUAgNTU\nVJM/aV5eHnx9fQ3barUaaWlpddokJiZi7969OHLkSINX+F20aBEAYF+qQOqVBPx1RQ2WJCj+M66I\nqCEcdRgtNTXV5K/lTb7KVlRUYNWqVThw4AAkScL999+P2bNnw8nJqcVP2til3m+ZO3culi1bBkmS\nIIRocNrqVnhs2ACk7gNyz5YBcG1xbURk4S5f/v3MYGqW6OhoREdHG7YXL17c6n02GR5Tp05Fhw4d\n8MILL0AIgc2bN2PKlCn4/PPPW/ykPj4+0Gg0hm2NRgO1Wl2rzbFjxxAXFwcAuHLlCnbt2gUHB4cG\nV3r5+en/zc3WtrguIlKACxeAO14vyPyaDI/MzEycPn3asD18+HCEhYW16kmjoqKQlZWFnJwcdOvW\nDVu3bsWWLbVP8Pvll18MX8+YMQPjxo1rdImwITzyOWVFZNUYHhahyZME+/fvj0OHDhm2Dx8+jLvv\nvrtVT2pvb4+VK1ciJiYGYWFheOKJJxAaGorVq1dj9erVLdrnrUMoF4rbQcvBB5H1unDh9//wJBtJ\nNLEGNiQkBOfOnYOvry8kSUJubi6Cg4Nhb28PSZJw8uRJc9Vax63jIbf08y3CXTevIemnAC4BJ7JW\nTz0FPPQQMHWq3JUo1p2vnS3R5BxPUlJSq57AnE6sOQK88w7gvkfuUoiorXDayiI0GR63zsW4fPly\nrY+h9bt1kMGS+PkBublyV0FEbYnhYRGaPOaxfft2BAUFISAgAA888AD8/f0xevRoc9RmPF9ffXi0\ncjhGRBZKCCAvj+FhAZoMj9deew2HDh1Cr169kJ2dje+++w4DBw40R23Gc3UFnJz0JxERkfUpKgLa\ntdPfSFZNhoeDgwM8PDyg0+mg1WoxbNgwHD161By1tQynroisl0bDUYeFaPKYh5ubG8rKynD//fdj\n8uTJ8PT0hIuLizlqM1pNDXDWbTCKk67j/pZft5GILBWPd1iMBpfqZmVloaCgAJGRkXBycoJOp8Om\nTZuQm5uLMWPGICoqyty11nHncrPCQv1H0nZyrkTJjZZfPoWILNRHHwE//AC08Hww0jPFUt0Gp63m\nzp2LDh06oH379rCzs4ODgwOmT5+O8ePHm+S6KG3BwwNwdqhGaYUTrl2TuxoiMjmOPCxGg+FRUFCA\nvn371rm/b9++yM7ObtOiWkqSAL8uFQD0U6NEZGUYHhajwfAoLS1t8JtuP9/D0vipdQB4zJzIKmk0\nvDSJhWgwPKKiovDxxx/XuX/NmjWtvrZVW/Lr6QCA4UFklTjysBgNrrZ67733MGHCBGzatMkQFseO\nHcPNmzfx1Vdfma1AY0Xc64yhn+2Hm8sgAA5yl0NEpiKEPjx8fOSuhNDEhRGFEEhJScGpU6cgSRLC\nw8MxfPhwc9bXqAZXDAQEAMnJQM+e5i+KiNpGcTHQowfQyJQ6NU+bXxhRkiQMHz7cogKjWbp3189b\nMTyIrAenrCxKk2eYKxLPMieyPgwPi2K94fHrr3JXQUSmxJVWFsU6wyMoCDh7Vu4qiMiUOPKwKFYZ\nHjnu/fH5ga7IzJS7EiIyGYaHRbHK8PgoJQSP5/4TX37ODzMnshoMD4tileERHqE/vyMz/brMlRCR\nyTA8LIp1hke4/t9TP8pbBxGZiBA8YG5hrDI8QkMBCTqcy2+Pqiq5qyGiVsvNBTp2BDp0kLsS+o1V\nhke7dkAPr+uo0dnh3Dm5qyGiVsvMBMLC5K6CbmOV4QEAceNuYLbbFjjxM6GIlC8z8/f5aLIIjV7b\nytI1en2WykrAzQ24ehVwdDRvYURkWjNmAPfdB8THy12JVWjTTxJUPCcn/TWuOG9FpHwceVgc6w0P\nQP/HxjMFiZRNpwPOnOExDwtj3eHRuzdw6pTcVRBRa+Tm6ldZdeokdyV0G+sOD448iJTv9GlOWVkg\nWcMjKSkJISEhCAoKQkJCQp3HN23ahH79+qFv374YPHgwTp48adT+D1VHYdG+Ydi/31QVE5HZ8XiH\nRZItPLRaLebMmYOkpCScPn0aW7ZswZkzZ2q16dGjB/bv34+TJ0/i9ddfx9NPP23UcySd6Y7Fxc9j\n13+rTVk6EZkTz/GwSI1+kmBbSk9PR2BgIPz9/QEAcXFxSExMRGhoqKHNoEGDDF8PHDgQFy5cqLOf\nRYsWGb6Ojo5GdHS0Ybt3PzsAwKm0GwA6mrR+IjKT06cBI984Um2pqalITU016T5lC4+8vDz43nad\nGrVajbS0tAbbr127FmPGjKlz/+3hcae+ffX/HjtpDyEASWpxuUQkB660Mok731gvXry41fuULTwk\nI17JU1JSsG7dOnz//fdGPUevXoCHSwUuXm2Pn3/Wf0YUESkIV1pZLNmOefj4+ECj0Ri2NRoN1PVc\nbvnkyZOIj4/H9u3b4ebmZtRzSBIwdJD+eAcPmhMpEA+WWyzZwiMqKgpZWVnIyclBVVUVtm7ditjY\n2FptcnNz8cgjj2Djxo0IDAxs0fPEv+SKD13+ghEheaYom4jMict0LZZs01b29vZYuXIlYmJioNVq\nMXPmTISGhmL16tUAgFmzZuGNN95ASUkJZs+eDQBwcHBAenq6Uc8zarQExPwCZKcAg58yeT+IqA1l\nZgJDhshdBdXDei+MeLsPPgBOngTWrGn7oojIdPr0AdatAwYMkLsSq8ILIzbX0KHAvn1yV0FExigs\n1B8wj4yUuxKqh22ER58+wJUrwMWLcldCRM21b59+yspettl1aoRthIdKBdx/P8Q+LrkiUoyUFGDY\nMLmroAbYRngAiC/4O3xmxnDwQaQUDA+LZjPhoVF1x8UbnXjog0gJLl3STzNHRMhdCTXAZsJj6BgX\nAMB339yQuRIiatK+ffqFLnZ2cldCDbCZ8Igdr/8j3PaVClVVMhdDRI3jlJXFs5nw6N0b6ON/DcU3\nnLBnj9zVEFGjUlKA2y7kR5bHZsIDACbNbA8HVOHsoSK5SyGihuTn65fW37osNlkkm1pA/cxzdpj1\nywK4u3QG8Fe5yyGi+qSkAA88oF9iTxbLpn47bm6A+x/HA5s3y10KETXkv/8FHnpI7iqoCbZxbavb\n6XRAQADwzTf6M8+JyHIUFwM9egDZ2fp3e9QmeG2rllCpgEmTOPogskSbNgFjxjA4FMD2wgMAJk8G\nNm4E1+wSWZh164A//lHuKqgZbDM8+vRBqnccxt+Tj/JyuYshIgDA8eP6aavhw+WuhJrBJsNDCOCv\nlX9D4gl/fLiiRu5yiAgA1q4FZszgKiuFsL0D5r/59lv9gg739pXIzndChw4mLo6Imq+yEvDxATIy\ngO7d5a7G6vGAeSs8+CAwJKIMxded8P471XKXQ2TbNm4E7r6bwaEgNjvyAPTnIg0fDnRyrkRWrhM8\nPExYHBE1z9WrQEgIsH07P27WTDjyaKVhw4ARA8vhVHkVOQcuyF0OkW16801g9GgGh8LY9MgDAC5c\nABw2fAyv7f8GDhwAHB1NUxwRNe3cOeC++4BTpwBvb7mrsRmmeO20+fAAoF9+FRurHzq//Xbr90dE\nzfOHP+jnjufNk7sSm8JpK1ORJGD9euCzz/TzrkTU9t5/H/j1V+D55+WuhFqA4XGLhwfwxRfAn/4E\n8eU2FBTIXRCRFfvySyAhQX+NOU4VKxLD43YDB0LsSsL/TM1HRPANnDwpd0FEVujAAeCZZ/TB4e8v\ndzXUQgyPO1T36Y+jfWbg0tV2uKd/NRLeqkYNT0InMo1t24AJE/TndURGyl0NtQIPmNejshJ47o8V\nWLfFGQAQFVKOd9e4YMgQkz8VkW24eROYPx/YsQPYupXLcmXGA+ZtxMkJWLvZGbt3A34e13H0Jxds\nnrJTP9xWbtYSmZ9WC2zZAkRE6D9eNiODwWElZAuPpKQkhISEICgoCAkJCfW2eeGFFxAUFIR+/frh\n+PHjZq5Qf+2rH8+3x9/+Wo2Xp5cA06cD99yjX877888AgKIi5glRHfn5wMqVQFgY8MEHwHvv6Rek\ndOokd2VkIrJMW2m1WgQHByM5ORk+Pj4YMGAAtmzZgtDQUEObnTt3YuXKldi5cyfS0tLw4osv4vDh\nw7WLb6Npq0YKB5KT9fO2iYlAp07odTEVpaIjIsNuIrifM4L73gV/f2DkSC4iIRty+TKQng6kpQF7\n9wJnzujP4Zg+XX8ehyTJXSHdxhSvnfYmqsUo6enpCAwMhP9vKy3i4uKQmJhYKzy2b9+OadOmAQAG\nDhyI0tJSFBQUwMvLS46S9ezsgJgY/W3VKlxPz0RpTEcUljljT5oz9qT93rRi7GOAn6d+CXDnzoC7\nO4SLK4YvGgqndhKcnSU4t5fg1E4FRycVVr1bBcnRAbC3199UKuiEhGXL9FeovnWTJP3DL75Ytzyt\nFvjXv/RfS9Lv/19VKmDWrPrb//vfde9XqYD4eNO0X7u2/vZ/+hPbW0x7IQCdDtqbNVi7VgDVNUDN\nb7ebN6FoL+tnAAAL40lEQVSqqcKfhv8ClJbqb4WFwKVL0OYXYG16H6DwMiCgv6hhwN1QDXgCf0rp\nVefdk8X014bb//ijfgbRFGQJj7y8PPj6+hq21Wo10tLSmmxz4cKFOuExbNgww9f+/v6GQDKXZ18G\nSkrcUFjYBVcKO+PapQ6oLnfA9nYquGZnwzkzE84VFXCuqIColJCadaLOPhxQhXc/c4dKp4OdTgeV\nTgdJCFTBEa/iZr3tZ/1PRwgAkCTcev9QBUfMqSitt/3kee76jdveAd4UDnjmenG97eNe7lzn/qqm\n2t/x7rJKOGBWeVG97Z94ue5VKK26/Uv1/zxnNfDzrK99dSPt4+a6w/DTFwISgCqdPWZVXq23/bRn\nXSDpdFAJAZVOB50k4aZ0F2bpKuptP2LzcFQ6OeGmkxOut2+PchcXlDp3xKxfX/m94U/6m90eLfLc\n3qyzH61Whb///fU699vZaZGXx/Zt2T4nJwe//PIrHB2HYvDgOt/SIrKEh9TMIeydw6r6vi8lJcUk\nNZne1jr31NQAqd8D5eVARYX+dvMmoNM5wunpytqNhYBdpRb/s7AGQgfotALaGgEhBFSSgNMb+fp3\njLd+RkLAsVrg2Vdv1LobAOxUgMvS87XvFAJONcDTr1yvU6edHeC69Fy99TfYfsnZetvH19deBbgu\n+8nG2p+tE641NUD8X2/UbW8n4Lqs/p9//Cv1tQdcEs7rN249hyTBUSshfsGtvysJ+nSRYGcvwWFF\nif4bf7upVCo41gDxz9bZPezsHBHwUXb99dTb3g4LFy6st/2lS2wvZ/vvvtNPnCxevLjuNxpJlmMe\nhw8fxqJFi5CUlAQAWLp0KVQqFRYsWGBo88wzzyA6OhpxcXEAgJCQEOzbt6/WyMPsxzyIiKyAYpfq\nRkVFISsrCzk5OaiqqsLWrVsRGxtbq01sbCw+/fRTAPqw6dSpk7zHO4iIyECWaSt7e3usXLkSMTEx\n0Gq1mDlzJkJDQ7F69WoAwKxZszBmzBjs3LkTgYGBaN++PdavXy9HqUREVA+eYU5EZGMUO21FRETK\nxvAgIiKjMTyIiMhoDA8iIjIaw4OIiIzG8CAiIqMxPIiIyGgMDyIiMhrDg4iIjMbwICIiozE8iIjI\naAwPIiIyGsODiIiMxvAgIiKjMTyIiMhoDA8iIjIaw4OIiIzG8CAiIqMxPIiIyGgMDyIiMhrDg4iI\njMbwICIiozE8iIjIaAwPIiIyGsODiIiMxvAgIiKjMTyIiMhoDA8LlpqaKncJbYr9UzZr7p81981U\nZAmP4uJijBw5Er169cJDDz2E0tLSOm00Gg2GDRuG8PBw9O7dG++//74MlcrL2v+A2T9ls+b+WXPf\nTEWW8Fi2bBlGjhyJc+fOYcSIEVi2bFmdNg4ODnj33XeRmZmJw4cP48MPP8SZM2dkqJaIiO4kS3hs\n374d06ZNAwBMmzYNX3/9dZ023t7eiIiIAAC4uLggNDQU+fn5Zq2TiIjqJwkhhLmf1M3NDSUlJQAA\nIQTc3d0N2/XJycnBAw88gMzMTLi4uBjulySpzWslIrJGrX3ptzdRHXWMHDkSly5dqnP/W2+9VWtb\nkqRGQ6C8vByPPvooVqxYUSs4gNZ3noiIWqbNwuPbb79t8DEvLy9cunQJ3t7euHjxIjw9PettV11d\njYkTJ+Kpp57C+PHj26pUIiIykizHPGJjY/HJJ58AAD755JN6g0EIgZkzZyIsLAxz5841d4lERNQI\nWY55FBcX4/HHH0dubi78/f3xn//8B506dUJ+fj7i4+OxY8cOHDhwAEOHDkXfvn0N01pLly7FqFGj\nzF0uERHdSSjUrl27RHBwsAgMDBTLli2Tu5xWy83NFdHR0SIsLEyEh4eLFStWCCGEKCoqEg8++KAI\nCgoSI0eOFCUlJTJX2jo1NTUiIiJCjB07VghhXf0rKSkREydOFCEhISI0NFQcPnzYavq3ZMkSERYW\nJnr37i0mTZokKisrFd23GTNmCE9PT9G7d2/DfY31Z8mSJSIwMFAEBweL3bt3y1GyUerr3/z580VI\nSIjo27evmDBhgigtLTU81pL+KTI8ampqRM+ePUV2draoqqoS/fr1E6dPn5a7rFa5ePGiOH78uBBC\niLKyMtGrVy9x+vRp8ec//1kkJCQIIYRYtmyZWLBggZxlttry5cvFk08+KcaNGyeEEFbVv6lTp4q1\na9cKIYSorq4WpaWlVtG/7OxsERAQICorK4UQQjz++ONiw4YNiu7b/v37RUZGRq0X14b6k5mZKfr1\n6yeqqqpEdna26Nmzp9BqtbLU3Vz19W/Pnj2GuhcsWNDq/ikyPA4ePChiYmIM20uXLhVLly6VsSLT\ne/jhh8W3334rgoODxaVLl4QQ+oAJDg6WubKW02g0YsSIEWLv3r2GkYe19K+0tFQEBATUud8a+ldU\nVCR69eoliouLRXV1tRg7dqzYs2eP4vuWnZ1d68W1of4sWbKk1uxGTEyMOHTokHmLbYE7+3e7bdu2\nicmTJwshWt4/RV7bKi8vD76+voZttVqNvLw8GSsyrZycHBw/fhwDBw5EQUEBvLy8AOhXqRUUFMhc\nXcu99NJLePvtt6FS/f5nZy39y87ORpcuXTBjxgz0798f8fHxuH79ulX0z93dHfPmzYOfnx+6deuG\nTp06YeTIkVbRt9s11J/8/Hyo1WpDO2t4vVm3bh3GjBkDoOX9U2R4WPPJgeXl5Zg4cSJWrFgBV1fX\nWo81dU6MJfvmm2/g6emJyMjIBs/PUXL/ampqkJGRgWeffRYZGRlo3759ncvuKLV/58+fx3vvvYec\nnBzk5+ejvLwcGzdurNVGqX1rSFP9UXJf33rrLTg6OuLJJ59ssE1z+qfI8PDx8YFGozFsazSaWsmp\nVLfOa5kyZYph+fKtc2IANHpOjKU7ePAgtm/fjoCAAEyaNAl79+7FlClTrKZ/arUaarUaAwYMAAA8\n+uijyMjIgLe3t+L7d/ToUdx3333o3Lkz7O3t8cgjj+DQoUNW0bfbNfS3eOfrzYULF+Dj4yNLja21\nYcMG7Ny5E5s2bTLc19L+KTI8oqKikJWVhZycHFRVVWHr1q2IjY2Vu6xWEQ2c19Kcc2KUYMmSJdBo\nNMjOzsZnn32G4cOH43//93+tpn/e3t7w9fXFuXPnAADJyckIDw/HuHHjFN+/kJAQHD58GBUVFRBC\nIDk5GWFhYVbRt9s19LcYGxuLzz77DFVVVcjOzkZWVhbuueceOUttkaSkJLz99ttITEyEk5OT4f4W\n988Ex2VksXPnTtGrVy/Rs2dPsWTJErnLabX/+7//E5IkiX79+omIiAgREREhdu3aJYqKisSIESMU\nuRyyIampqYbVVtbUvx9++EFERUXVWgppLf1LSEgwLNWdOnWqqKqqUnTf4uLiRNeuXYWDg4NQq9Vi\n3bp1jfbnrbfeEj179hTBwcEiKSlJxsqb587+rV27VgQGBgo/Pz/D68vs2bMN7VvSP1lOEiQiImVT\n5LQVERHJi+FBRERGY3gQEZHRGB5ERGQ0hgdRMxUVFSEyMhKRkZHo2rUr1Go1IiMj4erqijlz5shd\nHpFZcbUVUQssXrwYrq6uePnll+UuhUgWHHkQtdCt912pqakYN24cAGDRokWYNm0ahg4dCn9/f2zb\ntg3z589H3759MXr0aNTU1AAAjh07hujoaERFRWHUqFH1fmQzkSVjeBCZWHZ2NlJSUrB9+3Y89dRT\nGDlyJE6ePAlnZ2fs2LED1dXVeP755/Hll1/i6NGjmDFjBl599VW5yyYySpt9hjmRLZIkCaNHj4ad\nnR169+4NnU6HmJgYAECfPn2Qk5ODc+fOITMzEw8++CAAQKvVolu3bnKWTWQ0hgeRiTk6OgIAVCoV\nHBwcDPerVCrU1NRACIHw8HAcPHhQrhKJWo3TVkQm1Jz1J8HBwSgsLMThw4cB6K+mfPr06bYujcik\nGB5ELXTrMw9u/+yHOz8H4s7PRZAkCQ4ODvjiiy+wYMECREREIDIyEocOHTJf4UQmwKW6RERkNI48\niIjIaAwPIiIyGsODiIiMxvAgIiKjMTyIiMhoDA8iIjLa/wOiwbQ8oAVMYgAAAABJRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 5
    }
   ],
   "metadata": {}
  }
 ]
}