Coding problem

profile9qh0o1k52
dolo_rbc_linear_irf.ipynb

{ "cells": [ { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[33mwarning\u001b[0m: 29, 1 : Symbol e_z has no calibrated value.\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "\u001b[33mFutureWarning\u001b[0m:/Users/aaron/anaconda3/lib/python3.7/site-packages/dolo/algos/commands.py:34\n", " \n", "Panel is deprecated and will be removed in a future version.\n", "The recommended way to represent these types of 3-dimensional data are with a MultiIndex on a DataFrame, via the Panel.to_frame() method\n", "Alternatively, you can use the xarray package http://xarray.pydata.org/en/stable/.\n", "Pandas provides a `.to_xarray()` method to help automate this conversion.\n", "\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "Model:\n", "------\n", "name: \"Real Business Cycle\"\n", "type: \"dtcscc\"\n", "file: \"rby.yaml\n", "\n", "Equations:\n", "----------\n", "\n", "transition\n", " 1 : 0.0000 : z = (1-rho)*zbar + rho*z(-1) + e_z\n", " 2 : 0.0000 : k = (1-delta)*k(-1) + i(-1)\n", "\n", "arbitrage\n", " 1 : 0.0000 : 1 - beta*(c/c(1))**(sigma)*(1-delta+rk(1)) | 0 <= i <= inf\n", " 2 : \u001b[31m-1.8947\u001b[0m : chi*n**eta*c**sigma - w | 0 <= n <= inf\n", "\n", "\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAD8CAYAAABKKbKtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsnXl8m9WZ77+PJMt2vCW2kziJnZhshISwJRCgpQ3DsBQoW2mBUi4DtJQWhuksvcDcTofLHTq0vbdTpoUOlKWUDmVrgRQKacpSlhKamD2Qxdmd1bGNd0uW9Nw/XkmWZdmWZcuS7ef7+ZzPOe/Z3uc1wT+f7TmiqhiGYRjGSOPKtAGGYRjG+MQExjAMw0gLJjCGYRhGWjCBMQzDMNKCCYxhGIaRFkxgDMMwjLRgAmMYhmGkBRMYwzAMIy2YwBiGYRhpwZNpAzJJeXm5VldXD7ldZ2cn3d3dFBcXj7xRhmEYWU5NTc0hVZ06WL0JLTDV1dWsX79+yO2uuOIKfvWrX3HGGWdwxx13cOyxx6bBOsMwjOxERHYmU8+myFLg5z//OT/60Y9Yv349xx13HJdffjnbtm3LtFmGYRhZhQlMCuTl5fH3f//3bNu2jX/+53/mqaeeYtGiRfzd3/0d9fX1mTbPMAwjKzCBGQYlJSXcfvvt1NbWctVVV3HXXXcxd+5cbrvtNtra2jJtnmEYRkYxgRkBZs6cyT333MOHH37IGWecwb/+679SVVXFt7/9bXbs2JFp8wzDMDJCRgRGRM4SkU0iUisiNycozxWRx8Llb4lIdTj/dBGpEZEPwvFfxbR5Jdznu+EwbfS+yGHRokX85je/Ye3atZxxxhn8x3/8B/PmzeOiiy7iT3/6E3b3jmEYE4lRFxgRcQN3AZ8DFgOXicjiuGrXAE2qOh/4D+D74fxDwOdVdSlwJfBwXLvLVfWYcDiYto8YhBUrVvDYY4+xfft2brrpJl599VVWrlzJscceywMPPEBnZ2emTDMMwxg1MjGCOQGoVdVtquoHHgXOj6tzPvBQOP0kcJqIiKq+o6p7w/kbgDwRyR0Vq1OgqqqK733ve+zevZv77ruPUCjENddcQ1VVFd/61rd46623bFRjGMa4JRMCMwvYHfNcF85LWEdVA0AzUBZX5wvAO6rqi8l7MDw99i8iIoleLiLXish6EVk/Wju+8vPzueaaa3jvvfd4+eWXWblyJf/1X//FiSeeyNy5c7nlllt47733TGwMwxhXZEJgEv3ij//NOmAdEVmCM2329Zjyy8NTZ6eEwxWJXq6q96rqclVdPnXqoAdRRxQRYeXKlTz55JMcOHCAX/ziFyxatIgf/vCHHHPMMSxevJjbbruNTZs2japdhmEY6WBYAiMiL4rI2XF59w7SrA6oinmuBPb2V0dEPEAJ0Bh+rgSeAv6Hqm6NNFDVPeG4FXgEZyouaykpKeHKK6/k+eefZ9++ffzsZz9j+vTp3HrrrSxatIiFCxdyww038Lvf/c62PBuGMSYZ7gjmMOAmEfnXmLzlg7RZBywQkcNExAtcCqyKq7MKZxEf4GLgJVVVEZkMPAfcoqpvRCqLiEdEysPpHOBc4MNUP2q0mTp1Ktdddx2vvPIKu3fv5j//8z9ZsGABDz74IOeddx6lpaWsXLmSf//3f+ftt98mFApl2mTDMIxBkeHM+4vI2zgjhf/EGXF8BXhZVY8bpN3ZwI8BN/CAqt4uIrcB61V1lYjk4ewQOxZn5HKpqm4Tke8AtwBbYro7A2gHXgVywn3+EfgHVQ0OZMfy5cs1FV9ko4XP5+ONN95g9erVrF69mvfeew+A8vJyTjrpJE466SROPPFEjj/+eAoLCzNsrWEYEwURqVHVwQYTwxaYd1T12HD6b4B/BKaoamXKnY4i2S4w8ezfv581a9bw4osvsnbt2uhajcvlYunSpVHBOeGEE1iwYAEez4T2ZWoYRpoYLYH5uqreE/O8DLheVa9OudNRZKwJTDwNDQ385S9/4c0332Tt2rW89dZbtLS0AJCbm8vixYs56qijWLp0aTRUVFTQzwY7wzCMpBgVgRnrjHWBiScYDLJx40bWr1/PBx98EA379u2L1ikvL+eII45gwYIFzJ8/PxrmzZtn99sYhpEUyQqMzaGMI9xuN0uWLGHJkiW98g8dOtRLcDZt2hTdvRbLtGnTmD9/PtXV1VRVVVFZWUlVVVU0lJeX2+jHMIyksRHMOBrBDJW2tja2bt1KbW1tNGzZsoVdu3ZRV1dHd3d3r/p5eXlUVlZSUVFBRUUF06dP7xUqKiqYNm0aZWVlFBUVmRgZxjglq6fIROQs4E6cHV/3qeodceW5wC+BZUADcImq7giX3YLjqywI3Kiqq5PpMxETXWAGIhQKcfDgQXbv3s3u3bupq6uLpg8cOMD+/fs5cOAATU1NCdt7PB5KS0spKyvrFZeWllJcXExJSUk0xD4XFxdTWFjIpEmTcLnM2bdhZCNZKzBhZ5ebgdNxDlSuAy5T1Y9i6nwTOEpVrxORS4ELVfWSsFPMX+NsjZ6Jsx15YbjZgH0mwgRm+Pj9fg4ePMiBAweiobGxkYaGhj5xQ0MDTU1NtLe3J9V3QUEBBQUFFBYWRsOkSZOYNGkS+fn5CeO8vDxyc3N7xZF0JHi93oQhJycnGkzcDKN/snkNJursEkBEIs4uY8XgfODWcPpJ4Kdh32LnA4+G/Y9tF5Faek7sD9ankQa8Xi+VlZVUVia/Mz0QCNDS0kJLSwvNzc29QmtrK21tbbS3t9PW1hYN7e3ttLa20t7ezqFDh+jo6KCzs5OOjo5oeiRxuVy9BMfj8UTj2OB2u3ulI8+RdGyey+XC7XYnjOPTiYKIDBrHpxM9DxSAfp8HSycb91cWId3P8QxUPpy2yZQPxnDaD9b2xBNPZPr06Sn3nwyZEJhEzi5X9FdHVQMiEnF2OQtYG9c24ihzsD4Bx9klcC3A7NmzU/sCY1hEps9KS0tHrE9Vpauri66uLnw+X8K4q6uL7u5u/H5/v6G7uzthCAQCdHd3EwwGCQQCvUIwGIyWRYLP5+v1HAgECIVChEIhgsFgrzg+nSgEg0FUFVUlFAr1iQ1jqDz//POcddZZaX1HJgRmOM4u+8tPNJ+RcO5PVe8F7gUQkXoR2dm/qQNSjnM/TTZitqWG2ZYaZltqZNS2z33ucwMVD2bbnGTekQmBGYqzy7o4Z5cDtR2szz6oasrulEVkfTJzkJnAbEsNsy01zLbUmAi2ZWIlM2Vnl+H8S8W5UvkwYAHwlyT7NAzDMEaRUR/BhNdUbgBW0+PsckOss0vgfuDh8CJ+I45gEK73OM7ifQDHLU0QIFGfo/1thmEYRg8ZOcmvqr8Hfh+X992YdBfwxX7a3g7cnkyfaWawe28yidmWGmZbaphtqTHubZvQJ/nLy8u1urp6yO0++eQTOjo6mD59Om63e+QNMwzDyGJqamoOJbOGPaF9kVVXV5PKQcubbrqJH/zgB/h8Pm655Rauv/568vPz02ChYRhG9pHs7tsxc1xZRM4SkU0iUisiNycov05EPhCRd0Xk9fCp/7Tw/e9/n/Xr17N8+XK+/e1vs2DBAu69994+vrsMwzAmMmNCYMLuZe4CPgcsBi5LICCPqOpSVT0G+AHwo3TatGzZMlavXs3LL79MVVUVX//611myZAmPPvqoHXwzDMNgjAgMMe5lVNUPRFzBRFHVlpjHAvo5aDnSrFy5kj//+c8888wz5Obmctlll7Fs2TKefvppgsEBb2w2DMMY12SFwIjIcYNUSeReZlZ8JRG5XkS24oxgbuznXdeKyHoRWV9fX5+qyfF9ct555/Huu+/y8MMP09zczIUXXsi8efP4/ve/z6FD2XqQ2DAMI31khcAA3xikPBn3MqjqXao6D7gJ+E6ijlT1XlVdrqrLp05N+SB/QtxuN1/5ylfYvHkzv/nNb5g7dy4333wzlZWVXHXVVbz99tsj+j7DMIxsJisERlW/NkiVZNzLxPIocMFw7UoVj8fDRRddxEsvvcSHH37I1VdfzRNPPMGyZcs4+eSTeeSRR/D5fJkyzzAMY1QYdYERkVUi8mURKRhCs3XAUSKyLXy6/0biXMGIyPdE5CMReR94B0jVieWIsmTJEu6++2727NnDnXfeyaFDh7j88suZNm0aV155Jc899xx+vz/TZhqGYYw4mRjB/D/g08BHIvKEiFwsInmDtIlMhwk902UqIreJyHnh58PDcQjIA/aMpNHDpaSkhBtvvJGNGzeyZs0avvCFL7Bq1SrOPfdcpk+fztVXX80LL7xgW50Nwxg3ZOwkf3jr8V8BXwPOUtXiAeqeBNyqqmeGn28BUNV/76f+scBPVfVTA9mQ6Rst/X4/a9as4fHHH+fpp5+mpaWF0tJSLrroIs455xxOPfVUSkpKMmafYRhGIrL5RktEJB/4PHAJcBzw0CBNkrmkLJZrgOf7eXfWXDjm9Xo555xzOOecc+jq6mL16tU8/vjjPProo9x333243W5WrFjB6aefzhlnnMEJJ5yAxzOhnS8YhjGGGPURjIg8hiMOLwCPA6+o6oAnE0Xki8CZqvrV8PMVwAmq+rcJ6n4FuAH4bPhq5X7J9AimP/x+P2+++SZr1qxhzZo1rFu3DlWluLiYU089ldNPP52TTjqJpUuXkpOTk2lzDcOYYCQ7gsmEwJwFrIm42U+yTVJTZCLy18BPcMTl4GD9ZqvAxNPY2MhLL73EH/7wB9asWcOOHTsAyM/PZ/ny5axYsYITTzyRE088kVmz+hwPMgzDGFGyVmAARORkoJqYKTpV/eUA9T04U2SdOIv4BcBfx975IiLXAHfj3Adzqao+OZgdY0VgYlFVdu7cydq1a3nrrbdYu3Ytb7/9dnQn2qxZszj++ONZunQpS5cu5cgjj2TBggU2tWYYxoiRtQIjIg8D84B3gcgoRlU14cn7cBs3zrpLF47ATAJOw7mIbL2qrhKR14EjcHaZNQMfqOp5/XQJjE2BSYTP5+O9996Lik5NTQ1btmyJ+kTzer0cccQRUcE54ogjmD9/PnPnziUvb7ANfIZhGL3JZoH5GFisQ3jxUHaRicgvgGfH6wgmWTo7O9m4cSMffPABH374YTSuq6uL1hERKisrmT9/fjTMmzeP6upqKisrmTp1Ki5XVpzFNQwji8jmXWQfAhXAviG0GeousglPfn4+xx57LMcee2yv/KamJjZv3kxtbS1bt26ltraW2tpann76aeJ9s3m9XiorK6msrKSqqoqqqioqKyupqKhg+vTpTJ8+nYqKCgoLCxFJ5M3HMIyJTCYEphznkOVfgOgur0Gms5LyRZYM2bRNORNMmTKFFStWsGJFX31uaWmhtraWXbt2UVdXx+7du9m9ezd1dXW88cYb7NmzJ+FB0Pz8/KjgTJs2jbKyMsrKyigtLe0Tl5aWUlJSQlFRkY2ODGOckwmBuTWFNkP1RdYvqnov4fumly9fPnHvi05AcXExxx13HMcdl9i5dSgU4uDBgxw4cKBP2L9/PwcOHGDXrl28++67NDQ00NHRMeD7ioqKKCkp6RWKi4spLCzsNxQUFJCfn8+kSZOiIfbZ4/HYaMowsoRRFxhV/VMKzaK+yIjZRRZbQURygV/iXEr2KRFZr6o7hmmuEYPL5aKiooKKioqk6nd1ddHY2EhjYyMNDQ00NDTQ1NREc3NzwnDw4EG2bt1KW1tbNAz18jaXy0VeXl405Obm9krn5ubi9Xr7xJGQk5MTjfsLHo8Hj8fTJ+12u/F4PP3GLpcLt9vdJ8TmR9KxcXxaRExEjTHBqAmMiLSSeFpLcHaR9esqhgF8kRHeRQbcBpwTLq8A3gPMz0oGycvLY+bMmcycOTOl9qpKV1dXVGxaW1vp6Oigs7OTjo6OaIh97urqwufz0dXVlTD4/X5aW1s5dOgQfr8fv9+Pz+eLxt3d3dGQzTeTikhUbCLCExvi82Prx7eN5MWGoeYPFCL2DpQ3lOdE6cHKB2s3UD8jESebl0o/Q2kTm7722mtZtGgR6WTUBEZVi4bR/ATg/bhdZOer6ndj6hwDnK6qb4bPzewXERnKbjUjuxAR8vPzyc/PZ6Tv7kmGUCjUS3C6u7sJBAIEAoFoOjYOBoMEAgGCwWCvdKRNKBSKlkVCbF4kHRvHlqtqNC82qGrC8shzpCzyHJ+OzUtUbyhl8QFIKm+g8tjnROnBygdrl0zZcOJk81LpZyht4tOf+9znxo/ADJNkdpFF66hqQESagTKg3+ska2pqDolIqm79ywfqO8OYbalhtqWG2ZYaGbXtjDPOGKh4MNvmJPOOsSIwyewiS2qnWewuMuB/hRf9h26Qs8Yz6D7wTGC2pYbZlhpmW2pMBNuyap+oiJwlIptEpFZEbo4pqgOqROQzIvI28AugNK55CHhTRLaIyFU46y+N8e/QmCuTUxUXwzAMY3CyRmDEcQdzF84usMXAZSKyOFy8DliAIyJfxXEFsy6mbSkwA/gjznrNHcBrtv5iGIaRObJpiuwEoFZVtwGIyKPA+cBH4TWVG4AHcJxZ7gDqIrvIcLYtP40zalmHM132SprtzebRj9mWGmZbaphtqTHubUuLLzIRceHs+jpyCG0uxrnZMvbOlxWqekOCur8gxt+YiPwTkKeq/xZ+/hegU1X/b4K20TWYgoKCZansotjftp9mXzMzCmdQnDvQ7mrDMIzxR01NzSFVHXRrZ1pGMKoaEpH3RGS2qu5Kstlw3MEk3Tb+JH8qzi4ffOdBvvPyd9jSuoUTZp3Ad075DucuPNcOvxmGMSFIdvdtOtdgZgAbRORFEVkVCQPUH447mBFzJZMMVx17Fdtu3MY9595DfXs95z16HsfccwyPffgYwVDS96gZhmGMa9Lmrl9EPpsoX/txFRM+HLkZ556XPThrKV/WmEvFYur+gt5TZKVADRBxovU2sExV++wii2Uk3PUHQgF+/cGv+d7r32PjoY0sLFvIzZ+6ma8c9RVy3HadsWEY4w9J0l1/2kYwYSHZCBSFw8f9iUu4fgC4AVgNfAw8rqobROQ2ETkPQESOF5E64IvAPSKyIdy2Efg/OKK0DrhtMHEZKTwuD1ccfQUbvrmBJ774BJNyJnH1qquZ/ePZ3PzHm9ncsHk0zDAMw8g60jmC+RLwQ5zdXAKcAnxbk7gIbLRIx4VjqsoLtS/ws/U/4/dbfk9Qg5wy+xSuOfYaLl58MQXeghF9n2EYxmiT7AgmnQLzHo5vsIPh56nAH1X16LS8MAXSfaPlvtZ9PPTeQzzwzgNsadxCkbeIy468jGuOu4bjZx5vmwIMwxiTZIPAfKCqS2OeXcB7sXmZZrSuTFZVXtv1Gve/cz9PbHiCzkAn80vnc8HhF3D+ovM5qfIk3C532u0wDMMYCbJBYH4IHAX8Opx1Cc7ZmJvS8sIUGC2BiaW5q5nHNzzObzf+lhe3vUh3qJupk6Zy3uHnccGiCzjtsNPIz8kfVZsMwzCGQsYFJmzERcCncdZgXlXVp9L2shTIhMDE0uJr4fktz/PMpmd4bstztPhaKMgp4Mz5Z3L63NM5tfpUFpYttKk0wzCyimwRmAoct/ohYJ2q7k/by1Ig0wITiz/o55Udr/D0xqf53ebfUddSB0BFYQUrq1eycs5KTj3sVBaULjDBMQwjo2RcYETkq8B3gZdwRjCfxdk+/MAAbc4C7sTxN3afqt4RVx65FnkZ0ABcoqo7RKQaZ2vzpnDVtap63WA2ZpPAxKKqbG3ayis7XuHlHS/z8vaX2de2D4AZhTP4bPVnOWHmCSyfuZxjZxxLobcwwxYbhjGRyAaB2QScrKoN4ecy4M+qeng/9d04By1PxzmZvw64TFU/iqnzTeAoVb1ORC4FLlTVS8IC8+xQfJ9B9gpMPKrKlsYtvLLjFV7Z8Qqv7XotOsIRhCOmHsHymctZPmM5y2cu5+iKo5mUMynDVhuGMV5JVmDS6U25DmiNeW6l962U8fTrTTmmzvnAreH0k8BPZQLMF4kIC8sWsrBsIdcuc+5K29+2n5q9Nazfu56afTX8Yesf+OV7v3TqIxw25TAWT13M4vLFLJ66mCOmHsER5UdQlDucm6sNwzCSZ8QFRkT+IZzcA7wlIs/gOJ48H/jLAE2Hcy0ywGEi8g7QAnxHVV8b1odkORWFFZyz8BzOWXhONG9v617W713Pu/vf5eNDH/NR/Uf8Yesf8Af90TpVxVUcXn44cyfPZe4UJ8wrncfcKXOZnDc5E59iGMY4JR0jmMifyFvDIcIzg7QbzrXI+4DZqtogIsuAp0Vkiaq29HlJjLv+2bNnD2LS2GJm0UzOO/w8zjv8vGheIBRge9N2Pqr/iI/qP2JD/Qa2NG7htxt/y6GO3lduT8mbwtwpc5ldMpvK4spomFU0y4mLZ5HnyRvtzzIMY4wy4gKjqv87xabJeESO1KkLO8csARrDN1f6wu+vEZGtwEKcy8ji7evlrj9FW8cMHpeHBWULWFC2gPMXnd+rrMXXwvam7Wxr2hYNW5u2sqlhEy9tf4lmX3Of/sryy6gorGBawTSmF05nesF0Jx2OpxVMo2xSGWX5ZZTkleCSrLk01TCMUSZtazBh1zD/E1gCRP/sVdW/6qfJOmCBiByGM712KfDluDqrgCuBN4GLgZdUVcPvalTVoIjMxbleedtIfs94pDi3mKMrjuboisTee1p9rexp3UNdSx17Wpy4rqWOA+0HONh+kPV713Og7QCt/taE7V3iYkrelKjglOaXUppfSkluCZPzJlOSV9IrPTlvMsW5xRR5iyj0FlLoLTQPB4YxhknnIv9/A48B5wLX4QhDfX+VY65FXo2zTfmBiDdlYL2qrgLuBx4WkVqgEUeEAD4D3CYiASAIXDda3pTHM0W5RSzKXcSi8oFv/ezs7uRg+8Go8DR2NtLQ0eDEnT3x3ta9bKjfQHNXM82+ZkIaGtSGSTmTKPQWUuQtoii3iIKcAgq8BUzKmeSkc8Jpr5POz8kn35PfJ87z5JGf48R5njxy3blO7HFijyubbg83jPFBOrcp16jqMhF5X1WPCuf9SVUT3hOTCcbKNuXxiKrS5m+j2dfMJ12f0NzlxC2+Ftr8bbT6W2n1tUbjtu42Wn2ttHe30+5vp727nY7ujmi63d+OJn0Bal/c4ibXk4vX7SXXHY49ub3SXreXHFcOXrfXSbtj0q4cJ7gHjj0uDzluJ/a4PNG8+PxIcIu797PL3Sc/Ni8+bVOURjrIhm3K3eF4n4icg7OeUpnG9xljCBGhKNcZlVQWD/+fhariC/ro7O6kM9DZb+wL+OgKdOELhuO4Z3/Qjy/gwx8Kx0E/vqAvmt8Z6KTZ10x3sBt/0I8/6Kc75KS7g910h7qjcSAUGIGf1PAQBLfLHRWc/mKXuAat45ZwvQHSsf1E0tEyeur0KYvLGyh/sLJEdUay7lBDor4TBUHGnZeOdArMv4lICfCPwE+AYuBbaXyfMYERkej01xSmZNocwBG9QCgQFZ1IOhAKRJ9jy4MajOYlCsGQUx5bL5IXyQ+GgtHyROn+4pCGBiwPhnrXiX9fSEO96iTzHJuXKH84I9KxiiBDEi+R5Oon6vc/P/efnFx1clq/J20Co6rPhpPNwKkAImICY0wYRMSZGnPngN2ePWRUNSo4saE/gRpIqJKp2ys/LJz92RBpO1h5f2Xxtg21PL5MVQkxeNtYe3PduWn/b5hWZ5d9XiayS1Wz5vCJiNQDO1NsXg4cGrRWZjDbUsNsSw2zLTXGsm1zVHXqYJ2M9taZrJpgTOYH1B8isj6ZRa5MYLalhtmWGmZbakwE20Z7i8nEm1Q1DMOYoKTDF1kriYVEALuq0TAMY4Iw4iMYVS1S1eIEoUhVPeDc+yIim0SkVkRuju9DRHJF5LFw+Vthd/yIyAki8m44vCciF8a02SEiH4TLRuNwy72j8I5UMdtSw2xLDbMtNca9baO6yA/DvvdlEuAPn/qfAbwHzAw/7wCWq2rSi2bl5eVaXV09Yt82blCFQKD/EAo5IRjsHUeCau8wHEScEEm7XH3j+HR/we0eOB5nZxAMI13U1NQcysZFfhjGvS+q2hFTJ49hrulUV1cz4U7yh0Jw8CDs3Am7djkhko7EjQN42SkogOJiKCx00rFxYSHk54PXCzk5iUP8L/J4AQoGobvbEbLu7r5pvx+6unqCz9eT7ux0QkcHtLc7IRhM/meTl+d8W3+hpMQJkyf3H+fnm1AZ4x4RSWr3bSYEZjj3vhwSkRXAA8Ac4ApVjRyXVuAPIqLAPWGvyRMXVdi3DzZsgA8/7AkbNji/eGMpKoI5c2D2bDjxRJgxA8rL+4ayMshN/975EcXvdwSnowPa2npCa2tPHEm3tDjplpaesHMnNDc7oaVlcMHyeh2hiYQpU3qnY0Npae/n4mITJ2NckQmBGc69L6jqW8ASETkCeEhEnlfVLuBTqrpXRKYBa0Rko6q+2ufl4/U+mPp6eP11J6xb54hJU1NP+bRpsGQJXH01LFzYIyhz5jh/fY/XX2xeb88v/eGi6ghVczN88kn/cVOTE0fS27f35HV399+/291bfCICFEn3F6ZMAY856zSyj0z8q0z53pfYCqr6sYi0A0fieFveG84/KCJP4UzF9RGYcXEfjKrzS+u11xxBee012LTJKcvNhWXL4EtfgiOPdERlyRJHYIzhIeJMBxYUwMyZQ2+v6owem5r6hsbGnjiSrq93/rtGxGmg9aySkt6iU1bWN11W1js9ebIjaoaRJjIhMMO59+UwYHd42mwOcDiwQ0QKAJeqtobTZwC3jdL3jA4+H7z0Evz2t/D738PesCZPmQKf+pQzMvn0px1xGWvTWBMFkZ61qqqqwevHEgw6I6SIAMWHhobe6R07eoQq1M+1CCKOyMQLT6Ln2HRh4fgd8RojyqgLzDDvffk0cLOIdAMh4Juqeih8ydhTYU+kHuARVX1hdL8sDbS1wfPPw1NPwbPPOusDRUXwuc/Bqac6grJ4sbNwboxv3O6eEclQCIWc0U+sCDU09E03NDibPz7+2Em3Jr5EDnA2a/QnSgOJldc7vJ+BMeYY9W3K2URW3gfT0QFPPgm/+Q2sXu2MXMrL4YIL4MIL4bTTbIRipJ/u7oHFqL9ud0SFAAAgAElEQVR8v7//PouKhiZKZWXje31wDJMN98EYQ2HbNrj7brj/fucvzqoquO46uOgiZwrM5sqN0SQnB6ZPd0KyRNaYkhWk7dudeKD1pcjGh2RHSZF0vjkNyQZSEhgRmQ9MV9U34vJPAfaq6taRMG7cEwrBmjXw05/Cc885U11f+AJcfz2ccor95WaMLWLXmObMSb5dMOiITCJBihelXbvgnXecdGdn/33m5ydeQxoor7TUEVZjxEh1BPNj4J8T5HeGyz4/UGMROQu4E2cN5j5VvSOuPBf4JbAMaAAuUdUdInICPS4MBLhVVZ9Kps+sorkZHnoI7roLNm92dnh95zvw9a/DrFmZts4wRhe3u+eX/lDo7Ox/Gi8+/dFHPXmBAW4ajZ/G629HXmzaton3S6o/lWpVfT8+U1XXR/yG9UfYVcxdxLiKEZFVsa5igGuAJlWdH3YV833gEuBDHHcwUVcxIvI7nDMyg/WZebq64Cc/gdtvd0TmxBPhV7+Ciy+2dRXDGCr5+c4fZEP5o0zV2cAQL0D9iVMyu/Gg9zbxWBFK9BzJmwDClOrX5Q1QNtjkZzpcxSTTZ+ZQhcceg1tucf7Bnn023HorHH98pi0zjImFSI/rn6H4IQyFem8TjwhQU1NvkYqUxa4vDSRMxcWDH6RNFMbIH6SpCsw6Efmaqv48NlNErgFqBmk74q5iRCSZPjPD66/DP/4j/OUvcPTRzprLX/91pq0yDGMouFw9XhbmzUu+XbwwxW8Xjz1c29gIu3f3pAdySzRpUl9PD8mkR9kdUaoC8y2ccyeX0yMoywEvcGG/rRxG3FVMkn06HY+Wq5jaWrjpJudg5MyZ8OCDcMUVthvMMCYSqQqTquP7Ll6AYkUq1gPEli09+T7f4PaUlsLPfw6f/ezwv3EAUhIYVT0AnCwip+K4agF4TlVfSqJ5OlzFJNNnpF16XcX4/fDd78KPfuQcLLvtNviHf3DcixiGYSSDSI/37qFeKRLZ/BDveig+PdRDuykwrBUmVX0ZeHmIzUbcVQzwSRJ9pp8tW+Cyy6CmxnHdcvvtUFEx6mYYhjGBSWXzQ5oYF65iABL1Oaof9stfwje/6Sy+Pf00nH/+qL7eMAwj2zBXMcN1FdPS4hyM/NWvnPnMX/0KKitHxkDDMIwsJFlXMeYlcTisWwfHHQePPOKstbz4oomLYRhGGBOYVAiF4Ic/hJNPdpwCvvoq/Mu/2A4xwzCMGMb3MdJ0cd11zha/L3zBiadMybRFhmEYWYcJTCp89avOxV7XXmsOKQ3DMPphQi/yi0g9sDPF5uXAoRE0ZyQx21LDbEsNsy01xrJtc1R16mCdTGiBGQ4isj6ZXRSZwGxLDbMtNcy21JgIttkiv2EYhpEWTGAMwzCMtDChp8jKy8u1eqh+foyMouqEUKgnHR8i9WLrx+bFpxM990eiPR2xeZH0QHkDxQOlbT+JkS3U1NQcSmYNZkLvIquurmbYJ/mNPgQCjoODTz7pCc3NTtzSAm1tztXt8XF7O3R0OPeyJQp+f6a/LDtwu517qiIh/nmgkJMzcN5A6chzbDxQOpXn2HeaoGYvIpLU5qgJLTBGcqg6zlfr6mDPHjh4EOrre+LYcOiQIxiDkZPjOJguLOwdT5ni+OrLy+sbcnOdEPvLyOvt+8vJ7e4/uFy9RwTxz4P9HPoLkRFVKJQ4BIM98UDpQKAnnSgvEOh5jqQHygsEnLPAwaATBwKOWEfyE5XHl0XSA92blQ4SiU9/opRs+XDaDhYGe69rAi5ImMAYBAKwa5dzhU1trZOOiEldnRO6uvq2y8uDqVN7wuGHQ3m5IxKTJzuexidP7h2Kix0h8XpH/zuN4REK9RageBFKJEr9lSV6Hiw/mfL29oHLE9kwWrhcAwtRqsKY6qjxlFPS7+zdBGYCcfAgvP8+bNzo3CwQEZTt253/2SLk5DievisrYflyuOACJ11Z6eRPm+aEggKbxphIuFzOHwbj6Y8D1YFFcSghWcEcqrB2dAyt74Euwozl+efhrLPS+/M1gRmHdHfDpk2OmLz3Xk/Yv7+nTmEhzJ8PRx0FF10ECxY4z/PmwYwZE3M4b0w8RHr+sh8vREaa/YlWJG/OnPTbkhGBEZGzgDtx7m65T1XviCvPBX4JLAMagEtUdYeInA7cgXM1sx/4duQWTRF5BZgBdIa7OUNVD47C52SchgZ44w14/XUnvP12z62pXi8sXgxnnglHH+0IypIlMH26jT4MYzySTSPNURcYEXEDdwGn41x1vE5EVqnqRzHVrgGaVHW+iFwKfB+4BMd1wedVda+IHIlzwVjstW2Xq+q43xa2fbsjJK+95sQff+zke71w/PFwww1wzDGOoCxaNL7+OjMMY+yQiRHMCUCtqm4DEJFHgfOBWIE5H7g1nH4S+KmIiKq+E1NnA5AnIrmq6ku/2Zmju9sRkmefhd/9zlk/AWcR/VOfgiuucBbsli93Ft4NwzCygUwIzCxgd8xzHbCivzrhK5abgTJ6O1/7AvBOnLg8KCJB4DfAv+kYPkXa2Ogswj37LLzwgnOGxOuFv/or+Nu/hZUrnakuWysxDCNbyYTAJJr5jxeCAeuIyBKcabMzYsovV9U9IlKEIzBX4Kzj9O5Y5FrgWoDZs2cPzfI009oKjz8ODz/sjFiCQWe31oUXwuc/D6ef7izOG4ZhjAVS+vtXRPq9F1hEPj9I8zqgKua5EtjbXx0R8QAlQGPMu58C/oeqbo00UNU94bgVeARnKq4Pqnqvqi5X1eVTpw7q6SDtqMKf/gRXXunsSf/qV53dXjffDGvXwr598MADjsiYuBiGMZZIdQTzooicqao7YjNF5GrgfwG/G6DtOmCBiBwG7AEuBb4cV2cVcCXwJnAx8JKqqohMBp4DblHVN2Le6wEmq+ohEckBzgX+mOK3jQq7dsFDD8EvfgHbtkFREXzlK3DVVbBihe3wMgxj7JOqwPw9sEZEzlbVLQAicguOUHx2oIbhNZUbcHaAuYEHVHWDiNwGrFfVVcD9wMMiUoszcrk03PwGYD7wLyLyL+G8M4B2YHVYXNw44vLzFL8tbajCq6/CHXfA6tXO86mnwv/+385ZlEmTMm2hYRjGyJGyN2UROQ24B7gA+CpwPHCuqjaNnHnpZfny5Toazi5VnYX62293zqtMnw5f/zr8zd/AYYel/fWGYRgjiojUJHMhWcqL/Kr6ooj8DfAK8GfgNFVN4LFq4hIKwVNPwfe+5xx+rKqCn/wErrnGcehoGIYxnklJYESkFWdXlwC5wGnAQRERQFW1eORMHHsEAvDoo/Dv/w4ffeS4YLn/fmeNJRtO1xqGYYwGKQmMqhaNtCHjhZdeguuvdxxKHnkk/PrX8MUvOq7iDcMwJhJ2TG+E2LcPvvxlOO0052Ks3/7WcTB56aUmLoZhTEwyIjAicpaIbBKRWhG5OUF5rog8Fi5/S0SqY8puCedvEpEzk+0zXQQCcOedzl0ov/kNfPe78OGHzrkVO2VvGMZEZkw5uxSRxThblpcAM4E/isjCcJvB+hxx3nwTvvENZ6Ry5pnOAv6CBel8o2EYxtghE39jR51dqqofiDi7jOV84KFw+kngtPAGgvOBR1XVp6rbgdpwf8n0OWIcOuScuD/5ZCf95JOO3zATF8MwjB4yITCJnF3O6q+OqgaAiLPL/tom0yfg+CITkfUisr6+vj6lD/jWt5xT+N/+trOY/4Uv2Ml7wzCMeDIhMMNxdjnU/L6ZI+CL7Hvfg3fegR/8wPyDGYZh9EcmvCkPxdllXZyzy4HaDtZnH2pqag6JyM4hWd9DOb2vD8gmzLbUMNtSw2xLjbFsW1IXLmdCYIbj7HIV8IiI/AhnkX8B8BecEcxgffZBVVN2pywi65NxlZAJzLbUMNtSw2xLjYlg26gLzHCcXYbrPY5z+2UAuF5VgwCJ+hztbzMMwzB6yMQIBlX9PfD7uLzvRs6y4IjEfar6xdg6IvIZnJss5wCXqurzMcVTcUYyIZypNMMwDCODpOxNeaQJn4/ZTMxZFuCy2LMs4QOXxcA/AatU9clwfimwHliOs7hfAywbzLNzeXm5VldXD9nWruYuAp0BCqYWIG7bPmYYxsSipqbmUDJLDBkZwfRD9CwLgIhEzrJEBSZywZmIhOLangmsUdXIrZdrgLOAXw/0wurqalJx1//C37/AWz9+i9y2XE644QRW/N0KCqYWDLkfwzCMsUiym6OyyZlJ0mdZRrjtkDnrP87ia+u/xty/nstr33uNH8/5MS986wVa6lrS9UrDMIwxRzYJTNJnWYbTdiQOWgLMXDaTLz35Jb654Zss+dIS/vLTv3Dn3DtZ9bVVNGxpSLlfwzCM8UI2CUwy52OG3XYkDlrGMvWIqVzwiwu4sfZGll27jPcffp+7Ft3FYxc9xubnNhMKxM/mGYZhTAyyaZHfg7PIfxrOWZZ1wJcTbTcWkV8Az8Yt8tcAx4WrvI2zyN840DvTcWVy2/421t65lnfuf4eO+g6KZhZx9JVHc8xVx1C2oGxE32UYhpEJkr0yOWsEBkBEzgZ+TM9Zlttjz8eIyPHAU8AUoAvYr6pLwm2vBv453NXtqvrgYO9Lh8BECPqDbH5uM+8+8C5bfr8FDSlzPjOHY64+hsUXL8ZbYFdbGoYxNhmTAjPapFNgYmnd28p7v3yPdx54h8YtjXiLvBx+3uEcfv7hzD9rPrlFuWm3wTAMY6QwgUmC0RKYCKrK7jd2886D77DpmU10NnTi9ro57LTDWHTBIhZ+fiFFM+w2asMwshsTmCQYbYGJJRQIsfvPu9n49EY2PbOJpm3OmdDKEys5/PzDmXv6XCqOqcDlzqZ9GIZhGBkWmLBfsP8e7CR9psmkwMSiqtRvqI+Kzd71zga43OJc5nxmDnNWzqF6ZbUJjmEYWUGyApOuk/wVONcWvw08AKzWiTxUGgQRYdqR05h25DQ+853P0LqvlZ1/2smOV3aw45UdbH52M+AIzuxTZjPns3OYdcIsZhw7g9xiW78xDCM7SdsUWfiK4zOAq3B8hD0O3K+qW9PywhTIlhHMYLTubWXHnxyx2fnKTho2hw9yCpQtLGPm8pnRUHFMBd5C26FmGEb6yPQIhvD9LfuB/Tiu9acAT4rIGlX9n4naiMhZwJ30eFO+I648F/glsAxoAC5R1R1hJ5gfA5vCVdeq6nUj/1WZoWhmEUsvW8rSy5YC0H6wnb01e9m7fi/7avax45UdfPDfHwAgLqFsYRlTl0xl6uKeULawDE9eNrmeMwxjvJOuNZgbcS4MOwTcBzytqt0i4gK2qOq8BG2S8ab8TeAoVb1ORC4FLlTVS8IC86yqHjkUO8fKCCYZ2va3RUXnwLsHqP+4nsbaRjTo/PcVlzBl7hRHbA4vY8rcKdFQMrsEt9ed4S8wDGOskOkRTDlwkar28ripqiERObefNoN6Uw4/3xpOPwn8NDwVN+EprChk4TkLWXjOwmhewBegcUsj9R/V94QN9dSuriXoC0briUsoriruJTjFlcXRUDSriNziXOxHbRjGUEiLwKjqdwco+7ifokQekVf0Vyd8M2YzEPG/cpiIvAO0AN9R1dcSvURErgWuBZg9e/YgXzK28eR6opsHYtGQ0rq3laZtTU7Y3sQn2z6haVsTm5/dTPuB9j59eQu9UbEprCikYFoBBdMLKJzeky6Y5gRPrk3FGYaRXffBJOMRub86+4DZqtogIsuAp0Vkiar28Z+vqvcC94IzRTZMm8ck4pLo6GTOZ+b0KQ/4ArTta6OlrqUn7Gmhta6VlroWdv95N+0H2unu6E7Yv7fQS35pPvll+UwqmxRN55flk1+aT15JHnmT88gtye2Ttqk6wxg/ZJPAJOMROVKnLuwcswRoDG+B9gGoao2IbAUW4txyaQwRT66HydWTmVw9ecB6/nY/7QfaaT/YTtuBtmi6s7GTzoZOOho66GzspHlXMx0NHXQ1daGhgTXdk+fBW+TFW+gltyi3TzqnIAdvgRPnTOpJewu85Exy8jz5HnLyY+I8D558D548j03zGcYokk0Csw5YICKH4XhTvhT4clydVTibB94ELgZeCu9Wm4ojNEERmQssALaNnukTE2+BF+9cL1PmTkmqvoYUX4uPruYuuj7pwtccl/6kC1+LD3+bH3+rH1+rD3+rPypS/lY//jY/3R3dBP3BwV+YAHeu2xGcPA+eXE9POs+DO9eN2+vGk5s4HY1zwnFMcOW4cOckF7s8CdIeV/Q5NohLTBSNMUvWCEx4TeUGYDU93pQ3xHpTBu4HHhaRWqARR4QAPgPcJiIBIAhcN5irfmP0EZeQN9mZEqPvzNyQCHYH6e7opru9m+6Obvztfifd2U2gM9BvHPAFCHQFCPqCBLoC0RB97gzQ9UkXQV+QoD9IwBfolQ51h1IWt1QRt/QWHnecCEXK3YOk3TH13QPHA5a5kshzSZ90bJ1oOr4s5nnAOoPk9clPUJ6oHoIJ+ghivsjGyTZlY/RQVUKBUFRsoqE76OQNFIfbhQK9n4PdQTQY7jcuROr1Kg86sQZ652lQe8oSpKN9xNcNae86/cQa6p0ejwxFkAarO2D9JNsn0x8ukntHTDjmymMoW5jaHVWZ3qZsGOMWEXGmyXLc5EzKybQ5GUNVo8IULzx90rECFt8m0XN8/Zg6idolahvJQ+lbPlCfCcr6fYcmrtdf+/i+etnWnaC9Ju4vtu1gffcX5nxmTsoCkywTegQjIvXAzkErJqYc5yBpNmK2pYbZlhpmW2qMZdvmqOqgd85PaIEZDiKyPpkhYiYw21LDbEsNsy01JoJt5vvdMAzDSAsmMIZhGEZaGDMCIyJnicgmEakVkZsTlH9GRN4WkYCIXDwKJt07Cu9IFbMtNcy21DDbUmPc2zYm1mCS9LRcDRQD/wSsUtUnB+u3vLxcq6urh2xP+8F2fM0+iquKzQW+YRgTjpqamkPJLPKPld+Og3paVtUd4bJQsp1WV1eTyjmY9fes5483/RH/Rj/Lvr6MlbeupGBqwZD7MQzDGIuISFK7b8fKFFkiT8uzUulIRK4VkfUisr6+vj4lY5Z/fTk31t7I8uuWU3NPDT9Z8BP+/H//TMAXSKk/wzCM8chYEZhkPC0nhareq6rLVXX51KmDjvD6ZVL5JM7+6dl84/1vUHVyFWu+vYa7F9/NR7/5iLEw7WgYhpFuxorAJONpOSNMXTyVy39/OZe/cDmefA9PXPwED618iLq1dZk2zTAMI6OMFYGJeloWES+Ok8tVGbapF/PPnM91717HOT87h/qP67n/pPt54NMP8PFvP3bcNhiGYUwwxsQuMgARORv4MT2elm+P9bQsIscDTwFTgC5gv6ouGajPdDm79LX6ePfBd1n747V8sv0TJh82mRO/dSLHXHUMuUW5I/4+wzCM0SRZZ5ejLjAiUh3Z8RWTd7yqrhtVQ0i/N+VQMMSmZzbx5o/eZPcbu8ktyeW4rx3Hir9dQcnskrS91zAMI51ks8C8DXxeVfeEnz8L/FRVl46qIYyuu/66t+pY+x9r+ehJZ2f1grMXsPTypRz++cMntEdewzDGHtksMMcDdwOfB44DvocjOLsHbJgGMnEfzCc7P2Hd3ev44L8/oHVPKzkFORxx4REcedmRzD19Lu4cu5PeMIzsJmsFBkBETgLuwVkrOUdVUzuQMkwyeeGYhpSdr+3kg0c+4KMnPqKrqYv8snyWfGkJR152JFUnV+Fyj5U9GIZhTCSyTmBE5Hf0PruyGNgHNAGo6nmjYkgM2XKjZdAfpHZ1LR8+8iEbn9lIoDNAfmk+c0+fy7wz5zH/zPkUzSzKtJmGYRhAdgrMZwcqV9U/jYohMWSLwMTib/Oz+bnNbH1hK7Uv1NK2vw2AaUdOY96Z85h35jzmnDLHfKAZhpExsk5g+jXAcWR5qar+92i/OxsFJhZV5eAHB6ldXcvWF7ay6/VdBP1BPHkeZiybQeWJlVSeVEnliZUUzyrOtLmGYUwQsk5gRKQYuB7Hh9gqYE34+dvAu6p6/qgYEkO2C0w8/nY/O17ZwfYXt1O3to59NfsI+oMAFFcWR8Vm1gmzmLZ0GnkleRm22DCM8Ug2CswzOOstbwKn4RyI9AJ/p6rvjooRcYw1gYkn4Auw/9391K2tY8/aPex+czfNO5uj5cVVxUxfOp1pS6cx7chpTFs6jfJF5XhybXrNMIzUyUaB+SBy1iU8LXYImK2qraNiQALGusAkonVfK/tq9nHww4Mc/OAgBz88SP3H9YS6HXc14hZK55VSOr+UKfOmUDq/Jz3lsCm4vbZN2jCMgUlWYEbzT9nuSEJVgyKyfSjiIiJnAXfiuIq5T1XviCvPBX4JLAMagEviPQZMBIpmFFF0bhELz10YzQt2B2nY3BAVnYZNDTRubWTnazvxt/qj9cQlFFcVM7l6MsWVxT2hqiddMLUAcSVybm0YhtGb0RSYo0WkJZwWID/8LICqar+r1OERz13E3GgpIqtib7QErgGaVHW+iFwKfB+4JB0fMtZw57iZtmQa05ZM6/UTUVU66jto3NpIY60TmmqbaN7dzO4/76alriU68ongynFROL2QgmkFTphe0Ds9tYD8snzyS52QV5JngmQYE5RRExhVHc7cy6A3Woafbw2nnwR+KiKimd4ml8WISFQcqk6q6lOuIaW9vp2WuhYn7Hbi9oPttB9op/1gOwc3HKT9YDtBX7Cfl0D+lBjBmZJHbnFu31DixN5CrxMKvOQU5OAtcJ5zCnLMy4FhjDHGympvohstV/RXR1UDItIMlOGs9RgpIC6hcHohhdMLmblsZr/1VBV/q5+2A220H2ynq6mLzsbOhKGrqYvmXc34Wnz4mn342/z99huPK8dFzqQccvJz8OR78OR5ouloXq4Hd64bd647mu4Ve93R4MpxOemcmOccNy6PC1eOy4k9MXkeF+IWXO7+0+IWxNU7LS5BxEZxxsRjrAhMMjdaJnXrpYhcC1wLMHv27OFbZiAi0ZFI2YKyIbUNBUP42/y9BMff7qe7vbt3ut2Pv81PoDNAd2c3gc5Ar7S/3U/HoQ4CvgBBX7An7goQ8AX6TPWNNlGhSRTcjgBFnhF6CVOvPAmnU4mhTxp66kRtjSvvlRffLi6v3+9Pk8AOOkExWHF8e02uLNPlQ26boM6595zL7E+n93fgWBGYZG60jNSpExEPUAI0xnekqvcC94Kziywt1hpJ43K7yCvJc87s9J2lGzE0pAT9QYLdQULdoWg66O/9HOoOEQo4IdgddNLdPc8aVELBkBMHQn3TIe2pE5sOKqrq5MWGYEw6phylb1pjnnVoMSRIQ0+d6A8qrjwmL75dn7x+f/gDFKkOX3yGK25xxQMJZ5++Mlk+xLbxdXIK0u/FPeMn+ZMhLBibcc7P7MG54fLLqrohps71wFJVvS68yH+Rqn5pkH7rgZ0pmlVO9k6/mW2pYbalhtmWGmPZtjmqOnWwTsbECCa8pnIDsJqeGy03xN5oCdwPPCwitTgjl0uT6HfQH1B/iMj6ZPaBZwKzLTXMttQw21JjItg2JgQGQFV/D/w+Lu+7Meku4IujbZdhGIaRGLtwxDAMw0gLJjCpc2+mDRgAsy01zLbUMNtSY9zbNiYW+dNFeXm5VldXD7ldy+4WOho6yC/Lp6C8AE/+mJlpNAzDGDY1NTWHMr7IPxz/YSJyC477lyBwo6quFpHDgcdiupgLfFdVfywitwJfAyLXL/9zeN2mX6qrq0nF2eXO13ay7q51bHxqI8GDQWatmMVxXzuOIy85Em+hd8j9GYZhjCVEJKndt2kbwYT9h20mxn8YcFms/zAR+SZwVMzW4gtV9RIRWQz8GsdFzEzgj8BCVQ3G9b8HWKGqO8MC06aq/zdZG4frTbnjUAfv/+p93v7529R/VI+30MuSS5ew7GvLmHn8TDu9bRjGuCQbvCmn7D8snP+oqvqA7eGtxyfg3CUT4TRgq6qmeo5l2Ewqn8SJ3zqRFX+3grq1dbx939t8+MiHvHPfO5QvKmfxlxaz5ItLmLpkqomNYRgTjnQu8ifyHzarvzqqGgAi/sOSaXspzignlhtE5H0ReUBEpiQySkSuFZH1IrK+vr4+UZUhIyJUnVTF+fefzz/u+0fOvedcCisKee3fXuNnS3/G3Yvv5uXvvszBDw8OfurZMAxjnJBOgRmO/7AB24qIFzgPeCKm/GfAPOAYYB/w/xIZpar3qupyVV0+dWrK5yz7Jbc4l2XXLuPKl6/kH/b+A2fffTaFMwp57fbeYrN3/V7H7YdhGMY4JZ1TZMPxHzZY288Bb6vqgUhGbFpEfg48OwLfMCwKpxdy/DeO5/hvHE/bgTY+/u3HfPTER7x2+2u8+n9eZdLUScw7Yx7zzpzHvDPmUTi9MNMmG4ZhjBjpFJh1wAIROQxnMf5S4MtxdVYBV+KsrVwMvKSqKiKrgEdE5Ec4i/wLgL/EtLuMuOkxEZmhqvvCjxcCH47w9wyLWLFpr29n6x+2svWFrWz9w1Y++O8PAKg4toL5Z81n3pnzqFxRiSfPtj8bhjF2Ses5GBE5G/gxPf7Dbo/1HyYiecDDwLGE/YfFbAr4X8DVQAD4lqo+H86fhLM+M1dVm2Pe9TDO9JgCO4CvxwhOQoa7i2wk0JCy/9391L5Qy9bVW9n9592EAiHcXjczj5/J7E/PpupTVVSdXMWkskkZtdUwDAOS30WWlMCEf3m/CrymqhtHwL6sIBsEJh5fi4/tL29n1+u72P3Gbvau3xu9y2Tq4qmO2HyqipnLZlK+qByXx5wxGIYxuoy0wPwV8GngFJzDje8Cr6rqncM1NJNko8DE093Zzd51e9n1+i5HdP68G1+zDwBPnofpR02n4rgKZhw7gxnHzWDakdNsas0wjLQyogIT7tANHA+cClwHdKrqokHajOhJ/nD+DqA1nB+IfKSIlOKc8q/GmSL7kqo2DWTfWBCYeFcq4XkAAAmrSURBVDSk1H9cz/539rPv7X3R2NfiiI7L46J8UTlTF0+l/Ihyyo9w0mULy/DkmvAYhjF8RnoE8yJQgLMY/xrwuqoeHKRNWk7yhwVmuaoeinvfD4BGVb1DRG4GpqjqTQPZOBYFJhEaUpq2N0XF5sD7Bzj08SGatjdFN3eLS5gydwrlR5RTdngZU+ZOYcrcKZTOK6VkdglurzuzH2EYxphhpE/yv48zyjgS5zDkJyLypqp2DtAm3Sf54zkfWBlOPwS8AgwoMOMFcQml80opnVfK4osXR/O7O7tp2NRA/cf1HPr4EIc+PkT9x/Vs/cNWgr5gr/bFVcVR0SmZU0JxZTHFlcWUVDlp87FmGMZQSUpgVPXvAUSkELgKeBCoAHIHaJboNP6K/uqEb62MPcm/Nq5t5CS/An8QEQXuUdWIW+npkV1jqrpPRKYl823jmZz8HCqOqaDimIpe+RpSWve10rStKRo+2fYJTdua2PzsZtoPtPfpK29yXlR0CmcUUlhRSMH0AgorCnvC9EJyS3LNLY5hGECSAhO+rvgUnFHMTuABnKmyAZslyBuJk/yfUtW9YQFZIyIbVfXVQWzpeaHItcC1ALNnz0622bhCXELxrGKKZxUz55Q5fcoDvgCte1ppqWuhpa6F5t3NTnq383zggwO0H2gnFAj1aevOdTOpfFKvkF+W3/NcNom8KXnkT8knb3JeNG1TdIYx/kh2iiwf+BFQE/YZlgxpOcmvqpH4oIg8hTN19ipwIHLYUkRmAAnXiMIjnnvBWYNJ8lsmFJ5cT3S6rD80pHQ2ddK2v61XaD/QTsehDjobOuk41MH+d/Y7z40DzaZCzqQc8qbkkTc5j9ziXHKLc8krycNb7I2mc4tz8RZ58RZ6yS3KxVvo7RNyCnJwuW3rtmFkA8lOkf0whb5H/CS/iBQALlVtDafPAG6L6+uOcPxMCjYbSSIuYVKZMyKZtmTw2chQIERnUyedDZ10fdJFZ1MnXU3h+JOuaNrX7MPX7KOrqYvmnc10NXfha/HR3d6dtG3uXDfeAkdsciblRNPeAi+efA85+Tl48j1OelJO9Dman+eEnPycaDqan+vBneuOpj15Hlw5LpsWNIwEpG3fanhN5QZgNT0n+TfEnuQH7gceDi/iN+KIEOF6j+NsCAgA14d3kE0Hngr/z+wBHlHVF8KvvAN4XESuAXYBX0zXtxlDx+VxUTC1gIKpBSm1DwVC+Fp9+Nv8+Fv9Ttzm78kL53d3dONv99Pd3k13Rzfd7eHnjm46Gzvp7uwm0BlwyjqdOrEbHlLFnevuEZ/Y2Ovulef29hNy3bhzeue5clw9zzkxeTkxZTlJ5MXFLo8JojE6TOgrk8fLNmVjeGhICXQ5ohPwBQh0Bgh0OaG7szuaDnQGCPgCBH1B59nn5Eeeg/5gtDzoDzr5vkA0HfQHe+r4g73yA74Aoe4Qwe4gGkz//5MujyM0fQQonOfy9BakaFl8eVxefLrfOonK40OC+n2C+/+3d3YhUpVhHP/9XZa+DETrIjIyYyEjagsKwYhaulgt+rgru/TSwCD6kCDqoosuym4i6FOoiIi6Em9E7VbJXHXFSiNvQrQoqW50Z+bp4n1nPHN2ZnZnds+8r+vzg8PMec7H/PfPnPPfec/MeeZeRyPyQF1kcmg45jiXBVqmMFR27WhqKQA06o1W2LSC6GK9VWvMNEJtpnu9rdblsVEr1WrhdduW1drnaxdqbbW29brVZ2Z/GWTYaERzhlLbOvMMsuI2s7Yf6VDvUeu4nz730arNY17Lqg9dDxjHyYxlI+EEsJRu+dOoXwqecmhZ3dqWzQqs8rbl/RSWd9pXudaoFbaZmb1NcVlx3/WL9c777vH6zWVWt+z6P23evZmxjWOVvsbSeQc7jpMtzdDs+cu5JY41Qhi1hVC3cKpfCrhWvRS0XcOt3mV/pXVWja2q/G++oq/BSPqD8LueQbgB+HPOtdLg2gbDtQ2GaxuMy1nbrWY2Z0vgKzpgFoKkH+ZzkSsFrm0wXNtguLbBuBK0+S/SHMdxnErwgHEcx3EqwQNmcD6ce5VkuLbBcG2D4doGY8lr82swjuM4TiX4JxjHcRynEjxgBkDSpKSfJZ2K3TOzQdJpScckTUlKeh8cSZ9KOidpulBbKWmPpJPxsfstm4ev7Q1Jv0fvpiRtSqTtFkn7JZ2QdFzStlhP7l0Pbcm9k3S1pIOSjkRtb8b6bZIORN++ljT07nk9tO2U9FvBt/FhaytoHJF0WNKuOL9g3zxg+iS2gn4f2AjcCTwbWzznxCNmNp7BVyB3ApOl2qvAXjMbA/bG+RTsZLY2gB3Ru3Ez2z1kTU1qwItmtg5YD2yN77EcvOumDdJ7dwGYMLN7gHFgUtJ64O2obQz4G9iSkTaAlwq+TSXQ1mQbcKIwv2DfPGD6p9UK2swuAs1W0E6J2Ajur1L5SUJLa+LjU0MVFemiLQvM7IyZ/Rif/0s46G8mA+96aEuOBf6Ls6NxMmCC0JId0vnWTVsWSFoNPAZ8HOfFIvjmAdM/nVpBZ3GARZotpQ/F7p250dbaGsittfXzko7GIbQkw3dFJK0B7gUOkJl3JW2QgXdxmGeK0HBwD/ArcL7QKDHZ8VrWZmZN396Kvu2QlOpmOu8BLwPNO5OuYhF884Dpn/m0gk7JBjO7jzCEt1XSQ6kFXUZ8ANxOGMI4A7yTUoyk5cC3wAtm9k9KLWU6aMvCOzOrm9k4oQvuA8C6TqsNV1V80ZI2SXcB24E7gPuBlcArw9Yl6XHgnJkdKpY7rNq3bx4w/TOfVtDJKLaUBpotpXPirEJLa9SjtXUKzOxsPAk0gI9I6J2kUcIJ/Esz+y6Ws/Cuk7acvIt6zgPfE64TrVBoyQ4ZHK8FbZNxyNHM7ALwGWl82wA8Iek0Ych/gvCJZsG+ecD0T6sVdPxWxTOEds3JkXSdpOubzwktpad7bzV0mq2tIbPW1s2Td+RpEnkXx78/AU6Y2buFRcm966YtB+8k3ShpRXx+DfAo4RrRfkJLdkjnWydtPxX+YRDhGsfQfTOz7Wa22szWEM5n+8zsORbDNzPzqc8J2AT8QhjffS21noKutcCROB1PrQ34ijBcMkP45LeFMLa7FzgZH1dmpO1z4BhwlHAyvymRtgcJwxFHgak4bcrBux7aknsH3A0cjhqmgddjfS1wEDgFfANclZG2fdG3aeALYHmK91xB58PArsXyzX/J7ziO41SCD5E5juM4leAB4ziO41SCB4zjOI5TCR4wjuM4TiV4wDiO4ziV4AHjOI7jVIIHjOM4jlMJHjCO4zhOJfwPR4tMGQ3tgicAAAAASUVORK5CYII=\n", "text/plain": [ "<Figure size 432x288 with 8 Axes>" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "from dolo import *\n", "model = yaml_import(\"rby.yaml\")\n", "dr = approximate_controls(model, order=1)\n", "#dr_global = time_iteration(model, pert_order=1, smolyak_order=3)\n", "s0 = model.calibration['states']\n", "s1 = s0.copy()\n", "s1[0] *= 1.04\n", "sim = simulate(model, dr, s0, n_exp=100, horizon=40 )\n", "irf = simulate(model, dr, s1, n_exp=0, horizon=40 )\n", "\n", "\n", "print(model)\n", "i_z = model.variables.index('z')\n", "i_i = model.variables.index('i')\n", "i_n = model.variables.index('n')\n", "i_c = model.variables.index('c')\n", "i_k = model.variables.index('k')\n", "i_r = model.variables.index('rk')\n", "\n", "\n", "\n", "#plt.plot([1,2,3])\n", "\n", "plt.figure\n", "\n", "\n", "\n", "\n", "#for i in range(100):\n", "plt.subplot(811)\n", "plt.plot(irf['z']/model.calibration['z']-1, color='black')\n", "plt.ylabel('z')\n", "\n", " #plt.hold(true)\n", "plt.subplot(812)\n", "plt.plot(irf['i']/model.calibration['i']-1, color='black')\n", " #plt.hold(true)\n", "plt.ylabel('Inv.')\n", "plt.subplot(813)\n", "plt.plot(irf['n']/model.calibration['n']-1, color='green')\n", "plt.ylabel('Labor')\n", " #plt.hold(true)\n", "plt.subplot(814)\n", "plt.plot(irf['c']/model.calibration['c']-1, color='red')\n", "plt.ylabel('C')\n", "plt.subplot(815)\n", "plt.plot(irf['k']/model.calibration['k']-1, color='blue')\n", "plt.ylabel('K')\n", "plt.subplot(816)\n", "plt.plot(irf['y']/model.calibration['y']-1, color='purple')\n", "plt.ylabel('y')\n", "plt.subplot(817)\n", "plt.plot(irf['rk']/model.calibration['rk']-1, color='purple')\n", "plt.ylabel('Rk')\n", "plt.subplot(818)\n", "plt.plot(irf['w']/model.calibration['w']-1, color='purple')\n", "plt.ylabel('w')\n", "\n", " \n", " \n", " #hold\n", "\n", "#plt.subplot(212)\n", "#plt.plot(sim.iloc[i, :, i_i], color='blue', alpha=0.1)\n", "plt.show( )\n", "\n" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.23387445725364966" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model.calibration['i']\n" ] }, { "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.1" } }, "nbformat": 4, "nbformat_minor": 2 }