diff --git a/hank.ipynb b/hank.ipynb index e08b57d..017fd8d 100644 --- a/hank.ipynb +++ b/hank.ipynb @@ -455,7 +455,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "{'L', 'Div', 'Tax', 'A', 'NS', 'r', 'C'}\n" + "{'Tax', 'L', 'r', 'C', 'Div', 'NS', 'A'}\n" ] } ], @@ -483,7 +483,7 @@ "output_type": "stream", "text": [ "dict_keys(['nkpc_res', 'asset_mkt', 'labor_mkt'])\n", - "dict_keys(['Y', 'w', 'pi'])\n" + "dict_keys(['pi', 'Y', 'w'])\n" ] } ], @@ -931,18 +931,24 @@ "On iteration 0\n", " max error for nkpc_res is 0.00E+00\n", " max error for asset_mkt is 1.41E-01\n", - " max error for labor_mkt is 2.68E-02\n", - "Cannot solve constrained household's problem: No convergence after 30 iterations!\n" + " max error for labor_mkt is 2.68E-02\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ - "/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:26: RuntimeWarning: invalid value encountered in log\n", - "/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:4: RuntimeWarning: invalid value encountered in log\n", + "C:\\Users\\Bence\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:26: RuntimeWarning: invalid value encountered in log\n", + "C:\\Users\\Bence\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:4: RuntimeWarning: invalid value encountered in log\n", " after removing the cwd from sys.path.\n" ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Cannot solve constrained household's problem: No convergence after 30 iterations!\n" + ] } ], "source": [ @@ -1194,7 +1200,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.8" + "version": "3.6.5" } }, "nbformat": 4, diff --git a/het_block.py b/het_block.py index 9d05f1f..cf67800 100644 --- a/het_block.py +++ b/het_block.py @@ -282,7 +282,7 @@ def jac(self, ss, T, shock_list, output_list=None, h=1E-4, save=False, use_saved # if we're supposed to use saved Jacobian, extract T-by-T submatrices for each (o,i) if use_saved: return utils.extract_nested_dict(savedA=self.saved['J'], - keys1=[o.upper() for o in output_list], keys2=shock_list, shape=(T,T)) + keys1=[o.upper() for o in output_list], keys2=shock_list, shape=(T, T)) # step 0: preliminary processing of steady state (ssin_dict, Pi, ssout_list, ss_for_hetinput, diff --git a/krusell_smith.ipynb b/krusell_smith.ipynb index ed87860..12bd000 100644 --- a/krusell_smith.ipynb +++ b/krusell_smith.ipynb @@ -180,7 +180,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -252,14 +252,14 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 28, "metadata": { "scrolled": false }, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEKCAYAAAARnO4WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXd0VdeZvp9z1XvvXUICJIrovXcM2IDBFWNwSbVTJvYkv0ySmWScrElmMjOZycQ2BiGwsQ2mN9N7VQXU6KAuoX6Lbt+/P45sUIzM9QVhyn7W8rIl73PuFkt899z9vd/7KkIIJBKJRPLoo/m2NyCRSCSS+4Ms+BKJRPKYIAu+RCKRPCbIgi+RSCSPCbLgSyQSyWOCLPgSiUTymCALvkQikTwmyIIvkUgkjwmy4EskEsljguu3vYFbCQ0NFYmJid/2NiQSieShIS8vr0EIEebI2geq4CcmJpKbm/ttb0MikUgeGhRFue7oWnmkI5FIJI8JsuBLJBLJY4Is+BKJRPKYIAu+RCKRPCbIgi+RSCSPCbLgSyQSyWOCLPgSiUTymCALvkQikXxLCLug/XwT2kMV9+X1HqjBK4lEInkcEBY7hoJ6tEcrsda34xLkge/IGBS37n0GlwVfIpFI7hM2nRn9yRp0J2qw6y24RfkQ9ExPvPuGorh2/4GLLPgSiUTSzVjqDeiOVqHPrwOrwLNXML5jYvBIDkBRlPu2D1nwJRKJpBsQQmC63IruSCXG883gqsFnUAS+o2JwC/f+VvYkC75EIpHcQ4TVjuHsDXRHqrDU6NH4uuE/OR6f4VG4+Lp/q3uTBV8ikUjuAXaDBd3pWnTHq7G3mXEN9yZofiremeHd3ox1FFnwJRKJ5C6wNrajPVqFIbcOYbHjkRqI3/xUPNKC7uv5vCPIgi+RSCTfECEE5uttaI9UYSxpBI2Cd/8wfMfE4h7l821vr0tkwZdIJBIHETZBe3EDuiNVmCu0KF6u+I2Pw3dENC7+zp3PCyEwtLbgExh0j3f7VWTBl0gkkjtgN1rR59ShO1aFrcWEa4gngU+m4D0oAo27i1P3tFmtXDh5lLztmzDqdSz9r/fQaJy7l6PIgi+RSCRdYG0xojtWjf50LcJkwz3Rn8DZKXj2DkbROHc+b9TrOLdvF/mfb0XX2EBQdCxD5zyNsItuN7uRBV8ikUj+DnOlFu2RKtrP3QDAq28YfqNjcI/zc/qerfW15O/YwrkDe7AY24nv048pr/6ApMxBKJr7o+KRBV8ikUhQjcyMpU1oj1ZivtqG4uGC76gYfEdF4xro6fR9qy+UkrdtExdPn0DRKPQaOZZBs+YSnph8D3fvGLLgSySSxxq72YYhrw7dsWqsDe24BHoQ8EQyPkMi0Hg6VyLtdhuXTp8gd9tGai6ex8PHhyFz5pE5fRZ+waH3+CdwHFnwJRLJY4mtzYzuRDX6UzXYDVbcYn0Jfq4XXn1CUVycO583txsoOrCH/J1baK2vIzAiiolLvkPG+Mm4e3rd45/gmyMLvkQieaww1+jRHa3CUFgPdoFnegh+Y2JwT/B3elCqreEGBZ9v5dy+XZgMemJ6pTNu0SukDB7W7cqbb4Is+BKJ5JFHCIHpYgvaI5WYLraguGnwGRqJ36gYXEOdf/Kuu3KJ3G0buXDyKEII0oaNYtCsp4jq0dOh681GK6XHamiq0TPhxV5O78NRZMGXSCSPLMJix1BYj/ZoFdY6Axo/d/ynJeI7LBKNt5tz97TbuZyfQ972jVSWFOHu5cWA6bMZOGMO/mHhDt1D22Tk7IFKSo5UYTbaiE4NxGqx4eomdfgSiUTyjbDpLR1BI9XYdRbcIn0IWpCGd/8wp4NGLCYjxYf2k79jE8011fiFhjFu0Sv0nTgND2/H7I5vlGsp2FPO5bx6BJAyMIzMyfFEJPo7tadviiz4EonkkcFyoyNoJK8erHY8ewapQSMpgU6fz+uamyjctZ0ze3Zg1GmJTEnliR+9TdqwUWhc7vxELuyC60WNFO4tp+pCC26eLvSdEEu/ibH4h9zfRq4s+BKJ5KFGCIHpSiu6I1UYy5rAVcFnQAS+o6Nxi3DeyOzG9avkbd9M2bGD2Gw2egwezqBZTxHTM92hNw+r2UbZyVrO7Kugpc6Ab5AHI+f3IH10NB5e307plQVfIpE8lAibnfazDWiPVmGp0qHxccVvUjy+I5wPGhFCcP1MPrnbN3H9bAGuHh70nTSdQTOfJDAyyqF7GNrMnDtUSdGhKow6C2Hxfkx5JZ2UgeG4uNz+OEmYzSju3R+O0q0FX1GUa4AWsAFWIcTg7nw9iUTy6GNvt6I/XYPuWDW2NjOuYV4EzuuBz4BwFCebnlaLhbKjB8nbvomGiuv4BAUz+tmX6DdlBl6+jtkpNFXrKdxXzoVTddhsdhL7hpI5OY7o1K6Pk9oLC2lcmY2lvJzE9Z91u3/+/XjCnyCEaLgPryORSB5hrE1G9Xw+txZhtuOREkDgvFQ804KcNjJr12k5u2cnBZ9vRd/STGh8ItO//xN6jRqLi+udVTxCCCrPN1O4p4Ly4kZc3DT0GhlF5qQ4AiNu38gVNhvafftoylpJe0EBGj8/gp5ZiLBYuv0pXx7pSCSSBxrT9TZ0RyppL24E5YugkRjco32dvmdLbQ15OzZRdHAvVpOJxP4Dmf6Dn5LQN9Ohp2yb1c6l3DoK91XQUKHDy8+NobOT6DMuBq8ujpPsej0tGzbStGoVlooK3GJjifh//4/A+fPQ+Nyf0JTuLvgC2K0oigDeE0K8//cLFEV5HXgdID4+vpu3I5FIHgaE/ZagkXItiqcrfuNi1aCRAA+n71t1vpTcrRu4lHsSFxcXeo0ez+AnniI0PtGh6416C8VHqjh3oBJ9q5mgKB8mLOpF2tCILjX0lro6mj/8kOZP12Jva8MrM5Pwn/0Mv8mTUBxQ+dxLFCFE991cUaKFENWKooQDe4A3hBCHu1o/ePBgkZub2237kUgkDzZ2kw1Dbi3aY9XYmoy4BHviNzpGDRrxcK44/r2RmaePL/2nziRz2ix8g4IdukfrjXbO7K+g9HgNVpON2F5BZE6JJz49uMtPBMbSUhqzsmjbsRPsdvymTCH45cV4DxjQad0Nww3KmsoYEzvGqZ9PUZQ8R/uj3fqEL4So7vh3vaIoG4GhQJcFXyKRPJ7Y2szojlejO1WDaLfinuBP4MwkPNNDnD6fNxvbVSOzHZtvGpkt/S59xk3GzdMxu+Oay60U7i3nauENFI1C2pAI+k+OIzT29o1cYbejO3yYpqyVGE6dQvH2Juj55whetAj3uLhOa883nWdVySp2XN2Bl6sXBxYewMPF+U8vjtBtBV9RFB9AI4TQdvz3VOC33fV6Eonk4cNSq0d75KaRmVdGCL5jYvFIcH7yVNfUSMHnWzmzdycmvZ7otN7fyMjMbhdcKbhB4d5y6q624eHtyoBpCfQbH4tP4O0Lst1opHXLFppWZmO+cgXXiAjC3/oZgQsW4OJ/82cRQnC8+jjZxdmcqDmBl6sXC9IWsKj3om4v9tC9T/gRwMaOjzuuwBohxOfd+HoSieQhQAiB6XIL2sNVmC403zMjsxvXr5K7bSNlxw4j7HZSh45g0Ky5RKc5ZkpmNlopPV7D2f0VtDUY8Q/1ZMwzafQaEYl7F7741sZGmtd8TPPHH2NrasIjvTfRf/oj/tOno7jdVPmYbWa2X9nOqpJVXGq5RJhXGD8a+CMWpC0gwCPA6Z/5m9JtBV8IcQXo3133l0gkDxfCZsdwtgHd4UosNXo0vm74T03AZ1gULj5OGpkJwbUz+eRu20j5uULcPDzpP3UGA2c8SWBEpEP30DWrRmbFR6oxt1uJSglg1PxUEvuHouniOMl06RJN2dm0bt6CMJvxnTCB4JdfxnvokE5n+i3GFj49/ykfl31Mo7GRtKA03hn9DjMSZ+Dm4gZWExR8BDfKYOrvnPoz+CZIWaZEIulW7EYr+tO16I5VYWs14xruRdD8VLwzw1HcnDMyu+2g1HOL6T95Bp6+jsk1b1RoKdxbzqWceoQQJA8IJ3NKHJFJt3/iFkJgOHmSxqws9IePoHh4EDB3LsGLF+ORnNRp7bXWa6wuWc2Wy1sw2oyMihnF4vTFDI8arr4hGJogdzmcXga6OojoCxYjuDkfpegIsuBLJJJuwdpiRHe0Gn1OLcJkwyM5gMC5dzkopW3jzJ6dFO7a5tyglF1wvbiRwr0VVJ1vxs3DhT7jY+g/MQ7/Lo6ThNlM644dNK3MxlRWhktICKFvvkHQs8/iGnxT5SOEIK8uj+ySbA5VHMJV48rslNks6r2IHkE91EWNl+Hk/6lP9dZ2SJkEc9+F5AnQzVO2IAu+RCK5x5irdGgPV9J+7gYAXv3C8BsTi3uM84NSzbXV5O/Y7PSglNVi43yHkVlzrQGfQA9GzEshY3Q0Hl344ttaW2leu5bm1R9ira/HI7UHUe/8K/6zZqHxuNlgtdgt7Lm2h1UlqyhuLCbQI5DX+73Os72eJdQrFISA68fh+P/C+R3g4gb9FsLwH2AO6cW2s9Vc+Pw8P58hA1AkEslDgLALjBea0R2uxHSlFcXDBd9RMfiOisY10LljCiEE1edLyd220elBKUObmaLDVRQdqqRdayE0zpfJS9LpMbhrIzNzZSVN2atoWb8eYTDgM3IEUe/8Kz6jR3d6c9GatWy4uIEPSz+kVl9Lon8ivxr+K2anzMbL1QtsViharxb66nzwCoaxP4Mhr9HqGsyaU+WsPL6fujYTvSL9+PHkVDxlAIpEInlQ+TJR6kgl1vp2XALcCZiZhM/QSDRdKFvuhN1m4+LpE+Rt20jNJXVQathTC77RoFRzrZ7CvRWcP1mLzWonoW8ImZPjiUn7GiOzM2dozFqJdvdu0GgIeGImwS+/jGfv3p3WVeuq+bD0QzZc3IDeomdwxGB+OeyXjI0di0bRgLENTi+HU+9CawUEp8ATf4b+z1Ghg+UHrrI2twCD2caoHiH82/x+jEsL63bjNJAFXyKROMFXEqWifAh+pide/UJRunhyvhN3OyglhKDqQguFe8u5fk41Mus5IpLMSXEERd7eq0bYbOgOHKAxayXteXlo/PwIWbqEoBdfxC2ys8qnqKGI7OJs9lzfA8DUxKkszlhMRkiGuqClQi3y+avA1AYJo2DGHyFtOgWVrSxbV8rnRbVoFIU5/aN5ZUwSGdH3T5IJsuBLJJJvgLWhHe3RKgx5dQjLF4lSsXikBDj9hKptaqDg822cdXJQymazcym3nsK95V8amQ2ZlUTfcTF4+XVhZNbeTuumTeqg1PXruEVHE/GLnxMw/2lcfG++OdjsNg5WHmRV8Sry6/PxdfNlUfoinu/1PFG+Hf74Vflw4n+heJP6dcZcGPEDbFED2FNSxwfvnST3ejP+nq68PjaFl0cmEhnQvWqcrpAFXyKR3BHT9TZ0hytpL2kEjYL3gHD8xsTcVaJU/bUr5G3bSNnxI04NSpkMFoqPVHP2QCX6FhNBkd6Mf6EnPYdF4up++zcKa0MDzWvW0LzmY2wtLXj27UvMf/4ZvylTUFxvlkODxcCWy1tYXbKacm050T7RvDX4LealzsPX3RfsdijboRb668fAwx+Gfw+GfReDdxSf5VWyYs1BrjUaiAv24jez01k4OA4fj2+35MqCL5FIbouwC4wljWgPV6qOlV6u+I2PUx0r/Z1PlLrbQam2hg4js2M1WEw2YnoGMf6FniRkdO27Y7p8maaVK9VBKYsF34kTCVnyMl6DBnX6ZHLDcIOPyz5m7YW1tJpa6Rvalz8N/BOT4yfjqnEFswFyPoAT/wdNlyEgDqb9HgYsot7sTvaJa3x0aj8tBgsD4gN5e3ovpqZH4OrkMde9RhZ8iUTSCbvZhiGvDu3RKmyNqmNl4JwUvAdHoOniyflOWC0WSo8eIG/bJhory50alKq90krh3gquFNSjKAo9hoSTOSmesPgujMyEwHDqNI1ZK9AfOqwOSs3rGJRK6jwodaH5AquKVSMzq93KhLgJLM5YzIDwAeobgrYOcpZBznJob4LogfD0Cuj9JOdvtLNs6xW2FFZjsduZmh7B62OTGZTgWIP5fiILvkQiAcCmNaM7UY3+ZA12gxX3eD8Cpifh9TVPznfibgel7HbBtTMNFO4tp+Zyq2pkNjWevuPj8A26vdmYsFho+/xzGrOyMJWU4hIcTOgbPyTouee+Mih1q5GZp4sn81LnsSh9EQn+CeqiuhI48Vc4txZsFug5E0b+EBE3nKOXG1mWnc/hCzfwcnPh2aFxLB2VRGLo/QkzcQZZ8CWSxxxLXYdjZYHqWOmZHoLf2LtzrGyurSZv+2aKD+7Fav7mg1IWs43zJ2oo3FtB6412/EI8Gb0wld4jo7o0MrNptbSsXUfT6tVYa2txT04m8ne/JWDOnE6DUmabmR1Xd5BdnM2llkuEeoXy5oA3WZC2gEDPQHVQ6vJ+VT9/eR+4esHAl2D49zEHJLHlTDUfbDxKWa2WMD8P3prWkxeGxRPo3f0h5HeLLPgSyWOIEALTlVZ0hysxnu9wrBwSie/oGNycdKy8Z4NShyo5d6gKo85CeIIf017rQ3JmKJouzsEt1dU0rVpNy7p12PV6vIcNI/Kff4Pv2LEompvXtJpaWXdhHR+VfkRDewM9Anvwu1G/Y2bSTNxd3MFqhsKP1UZsXRH4RsDEf4LBr9CCLx+dKif7+H7qtSZ6Rvjxp6f7MSczGg/X+5tadTfIgi+RPEYIm532cw1oj1RhqdKh8XHDf0oCPsOdd6y0221cyjlJ7tYNXyZKOTModWZfBWUna7FZ7CT2C2XAlHiienQt92w/V0RTVhZtu3YB4D9jBsFLXsYrI6PTukpt5ZeDUu3WdkZEjeCdUe8wInqEeu/2Fsj7K5x6D7Q1ENYbnvwr9F3A9VYrK/ZcZW1uJe0WG2NSQ/nTgv6MTQ29L4NS9xpZ8CWSxwC70Yo+pxbd0WpsrSZcw7wImpeK9wDnHSstJiPFB/eRt30TLXU1BEREMnHJd+gzforDg1I1l1sp3FPO1bMNuLg4MChlt6M7eIimrCwMOTlofHwIfuklghe9iFt0dKe1526cY2XxSvaW70WDhhlJM1icsZiewT3VBc3Xbw5KmXWQNA7m/C/0mEReeTPLPi5id0ktLhqFOf1jeHVMEr2jnD/mehCQBV8ieYSxtpjQHa9Cf+oWx8qnUvDsGex0I9bQ2kLBru0U7t6OUdtGZI80Zj+/mB5DRziVKOXp48bgmYn0HReLdxdyT7vRSOvmLTStXIn56lVco6IIf/ttAhc8jYvfTZWOXdg5VHGIlcUrya/Px8/Nj8UZi3m+1/NE+nTIPqvy1PP5kk2gaKDPfBjxQ2wRfdldXMuyvx0nv7yFAC83vjsuhcUjE4nw/3YGpe41suBLJI8g5moduiNVGM7cAARefcPwGxODexdZrI7QVF1F3vaNlBzaj9ViJmXwMAbPmktMrwzHGrEmG6XHazizr5y2BiMBYV6Mey6NniOicOtqUKqpSU2UWrMGW1MTnunpRP/pT/hPn9YpUcpoNX45KHWt7RpRPlG8PeRt5qXOw8fNRx2UOr8Tjv/PzUGpET+EYd9F7xnButwKVqw+SHmTgfhgb/5lTgZPD4r91gel7jWP1k8jkTzGCCEwXWpBe7gS08UWFHcXfEdE4Ts6Btcg559Qq8pKyN22gUu5p3BxdSV97EQGPfEUITFxd74Y0LeaOHewkqJDVZgMViKT/e+cKHXlqpootWkTwmTCd9w4gpcu/UqiVJOxiU/KPuGTsk9oNjWTHpLOH8f+kSkJU9RBKUs75Gap0srGi50GperM7qw8fo2PThbRZrQyMD6QX8zoxdSMSFyc/PTzoCMLvkTykCNsdtrPNqD9IjrQzw3/6Yn4DotC4+WkY6XdxuWcU+Rs20DNhTI8ff0YPnchmdNm4RMY5NA9mqr1FO4r5/ypWuw2QXJmGJmT44lK6TpRqj03l8aslej270dxdyfgyScJfnkxHikpndZebb36ZaKUyWZiXOw4FmcsZnDEYPUNQd+gTsSeXgaGBojqD/OXQ/qTlNa3s2zLFbaeqcZmF0zLiOTVMckMSnDs53qYkQVfInlIsZus6E/XqdGBLaab0YEDwlFc76IRe2g/eds30lLb0Yj9ho6V1RdaKOhwrHR105A+Kpr+E+MIjPC+/TVWK9rdu2lckYWxqAiXwEBCv/99gp5/DtfQ0E73zq/PZ2XxSg5VHMJN48bslNm8lP4SyYHJ6qKGS3Dyr1C4BqxGSJ0GI99AJIzi8KVGPliZz5GLDXi7u/DCsASWjkoiPuT2+3oUkQVfInnIsLWZ0R2vQneyFmG04p7kT+CTd9mIbWulcNc2Cndtp/2LRuxPvkEj1mbncv4NCvaUc6Nci5efG0NnJ9FnXAxevrdvxNp0elo+W0fzqtVYqqtxT0gg8p9/Q8CTT6LxujkLYLVb2Ve+j+zibM41nCPAI+A2iVIn1PP5LxKl+j8LI36IKagHmwurWb7pKOfrtIT7efD29J68MDSBgC6Srh5lZMGXSB4SLPUGtIcrv5yI9eoTit/YWNzj7l0jNnnQUIbMnudwI9ZstFJ6rIYz+yrQNhkJjLizY6Wltpam1atpWbsOu1aL1+BBRPzy/+E7YUKnQSmDxcDGSxtZXbKaKl0V8X7x/NOwf2JOjzk3E6WKN6qFvioPvIJg7Fsw9DWalUA+OnWd7BMHuKFVE6X+Y0F/ZvePxt3JTz+PArLgSyQPMEIIzFfb0B6uxFjW9OVErN+YGFxDnJuIBag6X0ru1vVqI9bFRW3EzprrcCNW12zi3MEKig5XY263Ep0ayJhn00js07XvjrG0lMasLNp27AS7Hb9pUwlZsgSvfv06rbthuMGasjV8ev5TtGYtA8IH8NbgtxgfNx4XjQuYdJD7rnp001IOwcnwxH9A/+e51iZYvvcq6/LyMFrsjE0L488Lkxjd4+EclLrXyIIvkTyACLugvbgB7eEqLBVaND6u+E+Ox2dE9F1NxF7OPUXO1o5GrI8vw55ayIDpjjdiG6t0FO4p50JOHcIuSB4QTuaUOCKTum7E6o8epXHFCgwnTqJ4exP0/HMEv/QS7rGxndZebL5IdnE2269ux2a3MTlhMi+lv0RmeKa6QFurTsPmLgdjK8QNg2m/R6TNIK+ijfc/KWFPaR1uGg1PZkbz6phkekY6/+nnUUQWfInkAeIr1sQhngQ+lYL3QOetiS1mEyWH1InY5ppqAsIjvvFEbGVZM4V7yikvacLVXUPG2Bj6T4wjIOz2nzKE2Uzr9h00rViB6eJFXMPDCf/ZPxC4cCEu/v6d7n2y5iTZJdkcqzqGl6sXC9IWsKj3IuL8Oz5t1JWo/jZn14KwQa9ZMPINrNGD2VVcx7J3T1FY0UKgtxs/GN+Dl0YmEO73aAxK3WtkwZdIHgBsOjO6EzXoT1Zj11txj7t7a+KvNGJTUpn145+TOnQEGhfHowML9pTTWKnDy9+dYU8m02dsDJ5dfMpQHSvX0rRqNda6OjxSU4n6wx8IeGImivvN5q3FbuHzq5+TXZzN+ebzhHiG8MaAN1iYtvAWx8oD6vn85X3g5g2Dl8Dw76HziWdtTgUr1hyksrmdxBBvfvdkBvMHxeLtLkva1yH/dCSSb5EvMmL1uXVgtePZO1htxCb6O33m3FxTRd72TRQf3HezETtrHjG9HWvEmtqtlByp5uyBCnTNanTghEW96Dk0EpcufHcsNTWqY+Xatapj5fDhRP3r7/AZPbrTa2rNWtZfWM+HpR9SZ6gjOSCZ3478LTOTZ+Lh4qE6Vp75RC30dUXgE/6lY2Wd1btjUGofbUYrgxOC+NWsdCb3jnhkB6XuNd1e8BVFcQFygSohxKzufj2J5GHAVN6REVt8S0bs2Fjcwp3XhKuN2A1fWhOrE7FzCYl1rBGrbTJydn8FxUersRhtxPQMZNzzXx8daCwro3HFCrURKwT+06cTvHTJVxwra3Q1fFj6Iesvrkdv0TMschi/HvFrRseMRqNoOhwr/6aamWlrIKyXamTWbyEXGs0s236FTYVV2OyC6X0ieW1MMgPiH/1BqXvN/XjC/xFQCjzcNnMSyV0i7AJjWZOaEXutDcXTFb9xcfiOdD4j9otGbO7WjVRfKHWqEXujXEvh3nIu5dYjgB6DwsmcHEd4FwEoQgj0x4/TtHwF+uPHUby9CX7heYIWvYR7bEyntSWNJawsXsnua7sBmJY4jcUZi0kPSVcXtJTDyXchP7vDsXIszPkfRMokTlxt4v3VZzh4/gaebhqeHxrP0tFJJIQ8uIlSDzrdWvAVRYkFngDeAX7ana8lkTyoCKsdQ0E92sOVWG+04xLoQcCsZHyGRKBx0pzr7xux/mERTHj5O/Sd4HgjtrykicI95VSWNePm4ULfCbH0mxiLfxdyT2Gx0LZzJ40rsjCVleESFkrYT35C0LPP4BJwU6VjF3aOVh0luzib07Wn8XHz4cXeL/JC7xeI8o1SF1Xlq43Y4k3q133mw8gfYg3vy46iWt7/6zGKqtoI9XXnH6ak8eLwBIJ8HvxEqQed7n7C/y/gbUBqoySPHXaDBd2pWnTHq7BrLbhF+RD8bE+8+oaidJHedCfURux2Cndto13bRkRyKrN+/I+kDh3pWCPWaudiTh0Fe8ppqtbjE+DOiLkpZIyJxqOLyVObTkfLus9oWrUKa00N7ikpRL3zr/jPno3mlkasyWZi+5XtZBdnc6X1ChHeEfxs8M+YlzoPP3e/DsfKzzscK4+Cux+M+D4M+y46z0g+zalgRfZBqlraSQ7z4Q/z+jJ3QAyebg9PotSDTrcVfEVRZgH1Qog8RVHGf82614HXAeLj47trOxLJfcPabER3tAp9Ti3CbMcjNRC/hbF49Ah0vhFbW03etk0UH9qH1WwieeAQhsye/40ascVHqji7rwJ9q5mQGB8mvdyb1MERuHQxeWqpq6N59WqaP12LXavFe8gQIn/z669EB7YYW1h7YS1rStfQaGykV3Av/jDmD0xLnIabxg0sRshbqTpWNlwA/1iY+g4MfIl6sztZtzRihyYG889zMpjUK7xLJ02J8yhCiO65saL8AVhvqeEaAAAgAElEQVQEWAFP1DP8DUKIF7u6ZvDgwSI3N7db9iORdDfmKh3aI5W0n70BKHj3D8N3TAzu0b5O37P20gVytqznwunjuLi40HvMRAbPcrwRq2s2cXZ/BUVHqrAYbcT2CmLAlHji0oO7fKMwXrhA04osWrdvB5tNnYhduhSvvn07ravUVrKqZBWbLm2i3drO6JjRLM5YzLDIYR2OlY2qY2XOMtDfgMh+MPJNyHiKCw1Glh2Wjdh7gaIoeUKIwY6s7bYnfCHEL4BfdGxoPPCzryv2EsnDiBAC08UOD/pLHR70I2NUD/pAD6fvea0wj5wt66koOYeHtw9Dn3yaAdNnO5wR+/cTsT0GhTNgagJh8bc/XRVCYDh1isblK9AfOYLi5UXQM88QvPgl3OM6v7kUNxSTVZzFnut70Cgankh6gsUZi0kNSu148cvq03zhGrC2Q+rUDsfK0Zy42sSyVQUc6GjEPjc0nldkI/a+IXX4EokTCJug/dwNtIe+8KB3v2sPepvVQtmxw+Ru3UBDxXV8Q0IZt+gV+k2ahrvXneWaQgiqL7ZQsKfDmrhjIjZzUhz+oV00Yq1W2j7fRdOKFRhLSnAJCSHsR28S+OyzuAYFdbr3kaojrCxeSU5tDr5uvryc8TIv9H6BcO9wdVFlLhz7byjdqjpW9lsII97AGpLGjqJalv31OOeqWgnxceenU9JYJBux951uO9JxBnmkI3nQsZttGHLr0B6pxNashoH7jY29Kw96c7uBs/t2kbdjM7rGBkLjEhgyZz49R47BxfXOvjlfZMQW7Cmn/lobXn5u9B0fS99xsXj63v56u15Py/r1NK3MVq2Jk5IIXvKyak3scfOTicVmYftVtRF7qeUSEd4RLEpfxPzU+fi6+6qN2Iu74NhfoPw4eATAkKVqdKB7KJ/mVLD86FW1ERvqw6tjkpk3UDZi7yUPxJGORPIoYdNb0J+oRne8GrvBinu8H4GzUvDs7bwHvb6lmfydWzizewcmg5649L5Mfe2HJGYOcqgRazXbKDtRQ8HeCtputKsZsc/3pNfwr7Emrq+n+cOPaP7kE+xtbXgNGkTEP/0S3/HjOzVitWYt6y6s46OSj6hvryc1KJXfj/4905Omq41YqwnyV6uKm4bzaiN22u/VRqzJjZXHrvHhyTO0Ga0MSQySjdgHBFnwJZKvwdp0i+LG0mF9MC4Wj8Tbu0M6QlN1JbnbNlJyaB82m420oSMZPGceUT16OnS9UWfh3KFKzh2spF1rITzRn5FzU0jKDOs6I/byZXUidstWhNWK35QphCxdgldmZqd1tfpaPir9iHUX1qkTsVHD+O2o3zIyeqT6JtTeArkrVNdKXS1E9IV5yyBjLhcbjLy/5QqbC6ux2O1Mz4jktbHJDJSN2AcGhwq+oijzgH8DwgGl4x8hhJDTs5JHEnO1Du3hWxQ3mWH4jYvFLcL55mL1hTJytqxXrQ9cXekzYQqDnniKoKiYO18MtDW0U7i3gtLj1VjNdhL7hjBgajxRXcg9hRAYcnJoWpGF7uBBFA8PAp6eT8jLL+OekNBp7YXmC2QXZ7Pjyg4EgqkJU3m5z8u3TMRWwMm/3ZyITZ4Ac/+GSBrPyavNvH9LI/bZoXGyEfuA4ugT/h+B2UKI0u7cjETybSKEwHSlFe2hSkwXmu+N4sZu50pBDjlbNlBVVoynj+83DgOvv95GwZ5yLufVo2gU0oZGkDklnpAu5J7CakW7dy+Ny1dgPHcOl6AgQn/4QzUjNvimykcIQU5tDiuKV3xpTfxMr2dYlL6IGN+ON6Hac+r5fNF69es+81Vr4vA+t23Evjg8gWDZiH1gcbTg18liL3lU+TJs5FAllkodGl83/KclqIobJ3NPbVYLpUcPkbt1A42V5fiFhjFh8Wv0mTgVd887J1UJIagoaSJ/dzlV55tx93Qhc0o8/SbE4Rt0+zcfu8FAy4aNNK1ciaWyEreE+C4zYvde30tWcRYljSUEewbzxoA3eKbnMwR4BHRYE+9XC/2VA+DuC8O+C8O/h94rSm3EfjERG+rD7+f2lY3YhwRHC36uoiifApsA0xffFEJs6JZdSST3AWGxo8+vQ3ekCmtDuxo2MrcHPgPDUZwsXiaDnrN7Pyd/x2Z0zU2ExScy84f/QNqIMbi43vmv25ce9LvLaazS4RPgzsh5PUgfE41HF3JPa2MjzR99RPNHa7C1tuKVmUn4P76N38SJKLfYLfx9RmyifyK/GfEbZqfMVq2JbRY1ZOT4X9Qne98ImPQbGLyEeosXK49f48OOidghiUH8ZrZqTSwbsQ8PjhZ8f8AATL3lewKQBV/y0GFvt6I7WYPuWBV2nQW3WF+CX+iFV0ao04obXVOjqrjZsxNzu4H4Pv2Z9r0fk9BvgMNh4CVHqzmzT/WgD472YdLi3qQO6dr6wFxeTmNWFq0bNiLMZnwnTSRk6VK8Bw7stK6xvfHLjNhWUyuZYZm8NeQtJsRNUK2JTVo4/YF6Rt9aAaE9v7QmvviFNXGBbMQ+CjhU8IUQS7p7IxJJd2NtNamKm1O1CLMNj7Qg/MbG4pES4LTHTWNlBbnbNlBy+ADCbidt+CiGzJlPRHIPh67Xt5o4e6CS4sNVmAxqGPi453uS0Cekyz21nyuicflytLt3o7i4EPDUkwQvWYpHclKndddar5Fdks2WS1uw2C1MiJvAkj5L/i4j9l1VdWNshfiRMPPfEalTOHm1hWUfnmV/WT2ebhqeGaI2YhNDZSP2YcZRlU4s8D/AKNQn+6PAj4QQld24N4nknmCp06M9XIWhsB6EwKtfmJoqdRceN1VlJeRsXc/l3FO4unvQb/I0Bj0xl8CISIeub67VU7innLJTtQibIHlAGAOmJBCR9DUe9EeP0fjBBxhOnULj50fIK68QtOhF3MLDO60trC9kZfFK9pfvx03jxpwec3gp/SWSAjreEG6cV49tzq5Vj3F6z4ZRP8IaNZCdRbUs+78TnK1UG7E/mZzGohGyEfuo4OiRThawBljQ8fWLHd+b0h2bkkjuBaZrquLGWNqE4qbBZ2gkfmNicQ12LuBa2O1czjtNzpb1atiInz8jnn6OzGmz8PZ3TJdfc0m1Prh6pgEXNw3pI6PpPzmOwC6SroTFQtvnn9P4wXJM58/jGhFB+NtvE7hwAS6+N9+w7MLOwYqDrCxeSUF9Af7u/rza91We7/08oV6haiP22jG10F/4HFw9YcAiGPEDDH4JaiP2IzUjVjZiH10cslZQFKVQCJF5p+/dLdJaQXK3fJkqdagS8/U2NN6u+IyIVlOlugjevhNWi4WSw/vJ3baR5upK/MMiGDzrKfqMdzBsxC64eraBgt3l1F5pxdPHjT7jY+g3PhYvv9s/OdsNBlo++4zGlSuxVtfg3iOFkKWvEDDriU5h4Cabia2Xt5JdnM21tmtE+0TzUsZLzO0xF283b7DbVG+b43+BqjzwCoahr8PQ12gQfqw6fo1VJ6/TYrAwOCGI18cmy0bsQ0Z3WCs0KIryIvBxx9fPAY3ObE4i6Q6E1Y6hsCNVql5NlQqcnYz3kEg0XdgM3AmjXseZPTsp2LkFfUsz4YkpPPHmW6QNH+1Q2IjVYuPCKTVspKXOgH+oJ2OeSaP3yCjcPG5//VcUN4MHEfmrX+E7blwn64NWUyufnv/0Sw/63sG9+ePYPzIlYQquGlcwG+D0MtW1svkqBCXCzH+HzBe41ib4YM8V1uXmYLbZmdI7gu+MS2ZQgmNOnJKHF0cL/lLgf4H/RD3DP97xPYnkW8VusqI/VYvuaBW2NjNukXefKqVraiRvx2bO7NmJxdhOQr8BzPjBPxDft7/jYSOHqzizrwJDm5mweD+mvppByoAwNF3s6baKm1dewXvAgE7rqnRVrC5ZzYaLG2i3tjMqehRL+ixhaOTQWzzol8Hp98HQCDGDYPI/Q+/ZnKnS8v66MnYW1eCq0TBvYAyvjkmmR7jzvQzJw4WjKp1yYE4370UicRibzozuWDW6EzUIoxWP5ACC5qfikRbktOKmqbqSnC0bKD2yH7vNTtqI0ariJinFoev1rSbO7Kug+HAVZqONuPRgJk+NJ7Zn13tqP1dE44rlaHfdqrhZgkdycqd1pY2lZBVnsfvabhQUZiTNYHHGYnoGd/jvNF1Rn+YLPlI96NOmw8g3EfEjOHSxgfc+yOHElUb8PF35zrgUloxMJNzfuV6G5OHlawu+oihvCyH+qCjK/6A+2XdCCPFmt+1MIrkN1mYj2sOVGHLrEFY7Xukh+I2Pwz3O+djk2ksXOL3lMy6ePoGrqxt9Jk5j8CzHFTctdQYK9pRTdrIGYROkDApn4B3CRvRHj9G4fDmGkyfR+PoS8spSghYt6qS4EUJwqvYUK86t4ETNiS/DwF9Mf5FIn469VebB8Q4Peo3rlx70lpA0tp6p5v1NRymr1RLp78kvZ/bm2aFx+Hk618uQPPzc6Qn/CzsF2UmVfKtY6vRoD1ZiOFMPKHgPCFfNzLpQt9wJIQTl585wevM6yovO4OHtw7CnFjBg+myHPW7qrrVRsOs6lwtv4OKqKm4yp8QREPZ1iptdNC5fjqmsDNfwcMLfeovAZxZ2UtzY7Db2lu9lRdEKShpLCPEM4UcDf8TCngvxd/dXFTcXdqthI9ePqh70o36khoG7h/LJ6XJWHD1AdauRtAhf/mNBf2b3j8bdSb9+yaPD1xZ8IcTWjv80CCHW3fr/FEVZcJtLJJJ7iqm8De3BSowljShuGnxHROM7JtZpMzO73cbFUyc4vXkd9Vcv4xMUzNgXltBv8gw8vB1LlaoobSJ/l+px4+HtyqBpCfSbGIe3/9cpbtarHjfV1binpBD1+9/fVnGz+dJmVhavpEJbQbxfPL8e8WvmpMy5aX1w5lO10NcXg3+MGgY+aDE3zO6sPH6V1SdUD/phScG8M7cv43uGOX3EJXn0cFSWmS+EGHin790tUpYpgVtyYg9WYLrSiuLliu/IeyGt3EfOlvW01NYQFBXN4NnzSR87EVc3B1KlbHYuF9wgf9d1GipUj5v+k+PJGBONu2cXHjdNTTR/+OFNxc2gQYS88gq+4zsrbtrMbXxa9ikfln5Ik7GJPiF9WNp3KRPjJuKicQGzHvJXqWf0rRUQ1kt9ou/zNFeazSw7cpX1+ZVYbKr1wetjZRj448Q9k2UqijIDmAnEKIryl1v+lz9gdX6LEslXEXZBe1GHa2WVDo2/OwFPJOEzNApNFzLGO2EyGDizZwf5O7egb24iIrkHs3/yc3oMHYFG44C08otUqT3ltDUYCYzwZsKiXvQcGomLW9eKm6aVK2lZvwFhMuE7eRIhS1/Be2BnxU2dvo7VJatZd2EdBquBUdGjWNpnKUMih9xU3Jx+T1XctDdD/AhVWpk6lYLKVt77+By7Smpxc9Hw9KBYXhuTTJK0PpB8DXc6w69GPb+fA+Td8n0t8JPu2pTk8UJY7RjyOzT0De24hnoRND/1rnJi/z4+ML5vJjN+8FPi+zgmrTTqLRQdruLs/gratRYikvwZ9XQqSf26NlhrLyqmcfkHd1TcXGm5QlZxFtuubEMIwbTEaSzps4Rewb3UBc3X1Kf5/NWq4qbnTBj1Y+yxQzl4oZ53l53i9NUm/D1d+cH4HiwemUiYn3NHXJLHizud4Z8BziiKsgY15aoXqlrnvBDCfB/2J3mEsZts6E/XoD1Shb3NjFvM3btWttTVkrt1A0UH92CzWkkbOpIhTz5NZEqqQ9frmk2c2VdO8ZFqLCYb8RkhDJwWT3Rq16lS+mPHaVz+AYYTtyhuXlyEW8RXPW6WFy3nYMVBPF08WZC2gJfSXyLWL1ZdUHsOjv4XFG8ERQP9noFRb2IOSmXLmWreX3+YC3U6ogM8+dWsdJ4ZEoevh0wplTiOo78tU4D3gMuohT9JUZTvCCF2dtvOJI8sNr0F3XE1EFy0qxp6v6fT8OiiqDpC/bUrnN78GRdOHEXjoiF97EQGz55PcLRj8YHNtXoKdpdz/lQtQkCPQeEMnBZPaGwX0kqbDe3u3TQsW4appLRLxY1d2DlSeYQVRSvIr88nwCOA7/b/Ls/1eo5gz2BVcXP1sFroL+9Tw0aGfw+Gfx+tRzifnK5g+dED1LYZ6RXpx38+059Z/aJxc3KoTPJ442jB/zMwQQhxCUBRlBRgOyALvsRhrC0mdEcq0Z/uCARPD8FvfCwe8c5FIwshqCwt4vTmz7hWmIebpxeDZj3FoJlP4hsc4tA9aq+0kr/rOlfPNuDqqiFjTAyZk+PwD719KpXdZKJ102Yaly/HUl6Oe1ISUf/6O/znzEFzi+LGYrew8+pOsoqyuNRyiSifKP5xyD8yL3XeTY+b4k2q4qY6H3zCYdKvYfBS6i1eZB2/xocni9AarYxIDuHfnu7H2NRQqbiR3BWOFvz6L4p9B1eA+m7Yj+QRxFJvQHuoEkOB+itzt4Hgwm7nUt4pcjZ/Rs3F83j5BzD62ZfoP2Umnr53tgkQQlBe3ET+rutUX2zBw9uVwTMS6TehazMzm05Hyyef0Jidje1GA559+xL+l//Gb9Kkr6RKrb+4nlUlq6jV19IjsAe/H/17pidNx03jBhYj5GbB8f+BpssQlASz/hP6P8/lFivvb7/CxoIqrHY7M/pG8Z2xyfSLDXTqz0ki+XscLfjFiqLsANainuEvAHIURZkHMupQcnvMFVq0BytoL2lEcdXgOzwK3zExuAY5N9L/RU5szubPaKquJCA8gkmvfJ+M8ZNwc79z09Jus3Mpr578XWp8oG+QB6Oe7kH66K+RVjY00LRqNc0ff4xdq8Vn5EhC/vQnvIcN6/S0/UWq1Cdln9BmbmNQxCB+NfxXjIkZo65rb4Hc5XDyXdDXQ1QmLFgJveeQV9HGex8Xsae0DncXNWzk1TFJJIRIxY3k3uJowfcE6oBxHV/fAIKB2cioQ8ktCCEwXW5Be7AS06UWFE9X/CbEqRp6X+dCNMzGds7t203u9o3oGhvUnNg336Kng66VFrON0mM1FO4tR9toJCjKgfjAigoaV6ygdf0GhMWC37RphLz6Kl59Mjqtq9BWkF2czaZLmzDbzEyIm8DSvkvpH9ZfXdBWDSf/D3JXglkLKRNh1I8RiWM4cOEGf3v/FDnXmgn0duONiaksHpFAiK9U3Ei6BxlxKLknfOFD33agAkuFFo2fOwEzk/AZGommi6fnO2HU6SjYtZX8nVsxatuITe/D1Nd+SGLmIMdcKw0Wzh2s4uwBVVoZmRzAmIWpJPbtWgVkPH+exmUf0LZjR4e08imCly7BI6lzfGBpYylZRVnsur4LjaJhTsocFmcsJjmgQ4J544LqcXPmUxA2yJirpkqF92Xb2Rre/YvqcRMT6MVvZquKG293qbiRdC+ORhwmAW8AibdeI4To0kFTURRP4DDg0XHNZ0KI39zNZiUPHsIuaD93A+2BCiy1BlyCPQmc2wOfgREoXQwm3Ql9SzN52zdRuHsHFmM7yQOHMPSphcT07O3Q9YY2M2f2lXPuUBUWo42EPiEMnJZAdOrtz8KFELTn5dGwbBn6Q4fReHsT/PLLBC9e3Ela+YWZWVZRFserj+Pj5sPi9MW8mP4i4d4d6ypOq4qb89vVVKlBL8OIH9DuG8/a3Arezz5IVUs7aRG+/Hmh6nEjFTeS+4WjjxSbgOXAVsDu4DUmYKIQQqcoihtwVFGUnUKIk07sU/KAIax2DAX1aA91DEuFexG0MA3v/uEoLs4pSVrr68jZuoGiA7uxW230HDmGoU8+TVhC0p0vBtoa2inYU07psRrsNrvqWjktgbAunDSF3Y7u4CEaly2jvaAAl+Bgwn78I4Keew6XgJuRhTa7jX3l+1hetPxrzMx2qYW+/Dh4BsLYt2HYd2hR/Fl94jpZx/fTpDczKCGIf5mTwcRe4TJVSnLfcbTgG4UQf7nzspsI1aRH1/GlW8c/dzbukTzQCIsNfU4d2kOV2FpNuEX7EPxCb7wyQpwelmqsrOD05nWUHj2IomjIGD+JIXPmExQZ7dj11ToKdpVzIacORYFewyMZMDWBwIivca3csYPGDz7AdPESbtHRRPzqnwicNw+N1005psVmYduVbawoWsG1tmu3NzMr/FiND6wvAf9YmPYHGPgSNUYXlh+4yprTuRjMNib2Cud741MYkihTpSTfHo4W/P9WFOU3wG7UJ3cAhBD5X3eRoiguqJYMPYC/CiFOObtRybeL3WRFf7IW7ZFK7DoL7gn+BM7rgeddBI7UXr7I6U3ruJhzAld3dwZMn83gWXPxCwl16Pq6q23kfX6Nq2cacHXX0G9CLJmT4/DtQgVkb2+nZf0GmlaswFJdjUdqKtF/+iP+06ej3GKg9oW0Mrs4mzpDHb2De/Pv4/6dyfGTb5qZ5ayA4/8LbZUQng5z34M+87nUaOL9rZfZWFCFXcDsflF8Z1wKvaOcmzWQSO4ljhb8vsAiYCI3j3REx9ddIoSwAZmKogQCGxVF6SOEKLp1jaIorwOvA8THx3+DrUvuB3aDOhWrPdYxFZsaiP+EONyTApwq9F8MS53auJbrZwvw8PZh+NyFDJgxB2//AMeuP99M3s7rX9oTD36iQ0PfhQrI1tpK85o1NK1aja25Ga+BA4n4dUdO7C0/Q6uplTVla1hTuoYWUwuDIwbzLyP/hZHRIzuklc1qTuzJv0F7k2pmNuvPkDqVwspW/rbmDLtLVGnl80PjeXVMMnHBzvn1SyTdgaP2yGVAv7vxz+n4hKAXQvx7V2ukPfKDg01rRnu0Cv2JGoTZhmfvYPwnxjudLCWE4GpBLqc2rqX6QineAYEMeuIp+k+Z6ZgPvV1w9WwDeTuvUX9di3eAO5mT4skY27WG3lJXR9PKbFo+/RS7wYDvuHGEvP4a3oMGdVpXb6hnVfGqL10rx8eO55W+r5AZnqkuaKuBk39VB6bMOkidBmN+iogbxpGLDfzt4GVOXGnE39OVxSMTWTwykVAprZTcJ+6ZPfItnAEC+QbTtYqihAEWIUSLoihewGTg3xy9XvLtYG0xoj1UiT6nDmx2vPqF4T8hDrdI54aA7HYbF04e4/Smddy4fhW/0DAmLv0ufSZMcWhYymazcymnjrxd5TTX6PEP9WTc8z3pNSISV7fba/BNV6/SuHw5rZu3gN2O/8yZhLz6Cp49e3Zad73tOllFWWy5vAW7sDM9aTpL+ywlLShNXdB4WbU+OPMx2K2QMQ9G/wRbeAY7ztXw7qajFFe3EeHvwS9n9ua5YfHSzEzyQOPob2cEUKYoSg6dz/C/Ltg8CsjuOMfXAGuFENuc3qmkW7E0tKM9WIEhv8P+YGA4fuPjcOvCU+ZO2KwWSg4fIGfLZzTXVBMUHcv07/+EXqPG4eJ65187q9lG6fEaCnaXo20yEhztw5RX0ukxMBxNFzJGY1kZDe+9h/bzXSju7gQtWEDw0iW4x8Z2WlfWVMYH5z5gz/U9uCquzEudx+KMxcT5xakLas7C0f+Ekk2gcYMBL8LINzH6xbM+v5L3Vx/keqOB5FAf/ji/H08OiMbD1Tm/fonkfuJowf/G+nkhxFlgwB0XSr5VLLV62g5U0H72Brho8BkWid+4WFwDnbM/sJiMnNu/m5ytG9A1NhCelMLsn/6CHkOGOxQ4Ymq3UnSokjP7vhiW8mfMs2kk9ulaBWTIL6DxvffQHTqExseHkFdfJXjxS7iG3mz+CiHIq8vjg6IPOFZ1DB83H17OeJlF6YsI9epYd/04HPkzXNoD7n4w8g0Y/n3a3EL48OR1Vhw9QIPORP/YAH7x4kCmpEfiIqWVkocIRydtD3X3RiT3F3OFlrYDFWpWrLsLvmNj8Rsdg0sX5mF3wqjXcWb3DvK2b6Jd20Zs7z5Me/0NEvoPdKi5a2gzc3Z/BecOVWFutxKXHsyg6Qlf70N//DiN776HIScHl8BAVUP//PO4+Pt3Wne48jAfnPuAwhuFBHsG8+aAN3mm1zM3NfTnP1ef6CtOgncITPwnGPIa9RZPVhy9xkcnz6A1WRmTGsr3xmUyIiVEulZKHkocnbTVclND746qqdcLIaTW7CHDdL0N7f5yjOebUbxc8Z8cj+/IaDTezmXFtmvbyNu+mYLPt2JuN5A0YDBDn1pAbK+MO18MaJuM6rDU0WqsVjspmWEMnJ5AeMLtf7WE3Y5u/34a3nsf47lzuEZEEPGLnxO4YAGaW5q/VruVXdd2sbxoORebLxLlE8Uvhv6Cualz8XL1ApsVzq5TC319MQTEwYw/woBFXGsTvLfzCuvzKr90rfzeuBT6xNxZRSSRPMg4+oTfSZqhKMpTwNBu2ZGkWzBdbaVtXzmmSy1ofFzxn56I74goNE42GfUtzeRu28iZ3TuwmIykDhvJsLnPEJGU4tD1LfUG8ndd5/yJWgDShkUwcFoCQV00h4XVStuOHTS8/z7mS5dxi4sj8rf/QsBTT3XyoTfZTGy+tJmsoiwqdZUkByTzzuh3mJE046Y9cc5ydViq+RqE9oSn3oW+T1NS187/fVbGjnM1uGo0PD04ltfHJJMoc2IljwhO/W0XQmxSFOXn93ozknuLEALTlVa0+8oxXWlF4+umGpoNj0Lj7lyTUfv/27vvuKqr/4Hjr8PeewkyxI0KuLdpjpylDdOyskxtl+31/bX7Vo7SHOXemqWZVu5ciRtcgIsNIlP25p7fH59rSoEg3wCB83w8eHD53M/n3sO5lzefez7nvN/pqRzfspGzu3dQWlJC61596T5mLE6e3lU6Pv1KLie3R3PpeBIGhga06+tO4BAvbBxvUXDk582kLV5McXy8tlhqxgxsht6DuOnib05RDhsubmBV2CpS81Pp4NSB17u+zgDPARgIAyjIghNz4fB8LT2xR2cY8hm0Hs7JuAzmrTrFH+eTsTI1YnI/Xyb1boaLTfWuYyjKnaqqQzr33/SjAdAFlSbhjiWlpPByBll7YimKztIyVycW0OAAACAASURBVI701TJXVjPQZyYncXzLT5zbuwudTodf37vpNvqhKpcQTInN5uS2aCJOpWBkbEDAQE8CB3thaVv+1Exdbi7XfthA+rJllKSkYObvj+u772DVvz/C4MYsnfSCdFaHrWb9hfVkF2XTo0kPvuj7Bd3cumnj7DkpcHQBHFsMhZng2x/6LEb69OXg5TTmLTrK0ah07C2MeW1wKx7v6YNtNYe3FOVOV9Uz/FE33S4BooH7/vXWKP8TKSUFF6+RvSeWothsDG1NsLuvOZZd3KqdufJaYgJHN/9I+MG9gKD9gEF0u+9BbF3cqnT81chMTm6LJvpsGiZmhnQe6k3AQM9bropNX72aaytXUZqZiUWPHrh/9SUWPXqUuVCalJvE8tDl/HTxJwpLCxnoNZBJHSbR3qm9tkNGrFZVKngVlBRA21HQZxq6Jh3ZGXaVefOCOJuQiZuNVhB8fDeVnlhp+FQ+/AZASklBeDpZf8RSHJ+DoZ2plqK4syuiggIflUmLj+XIph+4EHQQQyMjAgYPp+u9D1Qpz42UkiuXMjjxezTx569hamlE93ub0aF/U0wrOHsuSUkhfcUKrq1dp62KHTAAp6lTMA8MLLNffHY8S88tZfPlzeikjhG+I5jUfhK+dvo89KmXtKmVZzdoP/uPg94vU+zQgi2nrrDghwNcTs7B29GCL+7vwJhOHmoOvdJoVHVI5yvgUyAf2A4EAK9IKVfXYNuUSkidpCAsTQv0V3IxdDDD/oGWWHR0qXagT46O5OimH7h4LAhjE1M6jxxNl5FjsLSzr7w9UhIXls6JbdEkXs7E3MaEXve3uHX6g4QE0pYsIeOnjciSEmyGDcNxyuR/rIqNzIxkydkl/Bb5GwbCgDEtxvBk+ydpaq1fVJV4Bg7OhLBftDz0XSdDrxcosGjCjyfi+G6/loe+jZs1c8Z3ZHh7N4xUHnqlkanqZ9ghUso3hRBjgHi0mrZ7ARXw64DUSfLPpZL9RyzFV/MwcjTD/qFWWAQ6I6oZxBIvX+DIph+IPHkME3Nzuo9+iE7D76taQrO/5bmxsjel78Ot8OvdBKMKrhkURkaRtnAhmb/+CkJgN/o+HCdNwsTHp8x+59PPs+jMInbF7MLU0JRH2j7CE35P4Grpqu0QdxwOzoCL27XFUn2mQY/nyDayY83RWBYf1BZLdfSy4+P7tDz0ag690lhVNeBf/xw+HFgnpUxXfzS1T+ok+WdSyPojjpLkPIyczXF4uDXm/s7VLjoSfz6UIxvXE3MmBDNLK3o99Cgdh47CzMqq0mN1OklEcDInt0WTlqDluRkwoQ2te7hVWCu28NIlUhd8R9a2bQhTU+wfGY/jk09i3KRJmf1Op5xm0ZlF7I/fj5WxFU93eJoJfhNwMHPQFktFHYAD07Xv5vYw4D3oNpl0nSXLD0WxPCiYrAJtsdRz/TvSw9dBBXql0atqwN+qz5iZDzynT4xWUHPNUm6mlRFMJWtPDCXJ+Ri5WuAwvg3mt6jNesvHk5K40LMc2biOuLCzmNvY0veRiQQOGY6JeeWZK0tLdVw6lsTJ7TFkJOVh72bBoCf9aNnlFnluwsNJXfAd2Tt3YmBhgePTk3CYOBEjR8cy7Tp+9TgLzy7kaOJRbE1teSHwBca3HV+2stSBGRB/DKxcYcin0PlJrhYYsWhPJGuPxpJfXMo97Vx5rn8LAjzLL2uoKI1RldIjAwgh7IEsKWWpEMICsJFSXv03G6PSI5cldZL80FSydsdSkpSHkYsFNoO8MG//vwX6wz+tJT78HJb2DnQddT/+A4dibFb5nPPSYh3njyQSvCOGrNQCHJta0WWYD807OlfYnvyzZ0mdv4CcvXsxsLLC/rEJODz+OEb2N64JSCk5mHCQRWcWcSrlFE7mTkxsN5GHWj2EhbEF6HQQvkUbo796RlsV2/tlbVVsZinfH4hg48kESqXkvgB3nunfnFau1UvjrCj1TU2kRwZoC/gIIW4+ZuVttUypEiklBaFpZO2OpfhqrjZ0M7415h0qDqyVPd7fA/2AiVPxH3gPRiaV584pKS4l7M8rBO+IJTejENdmNvQd2wrvDhXnlMkLDiZ1/gJy//wTA1tbnF56EYcJE8rkudFJHXti97DozCLC08NpYtmE97q/x5iWY/QlBEvg9Hpt1k3qBXBoDvfNB/+xnE/JZ8HGcLaevoKRoQFjuzZlar/mquCIotxCVWfprAKaA6eAUv1miQr4/yopJQVh6WTtjqE4MRcjJ/0YfUAdBfqiUkL/vELIjhhyM4to0sKWgU+0pWmb8ssaSinJO3ac1PnzyTt6VCsK/tqr2I9/BEOrG+kJSnQlbIvaxpKzS4jIjMDbxpuPe33MSN+RGBsaQ0khnFiqFQXPiAHX9vDgUvAbzamEbOauPs3u8CQsTQyZ3NeXSX3UqlhFqYqqnuF3AfxkVcd/lNsipaTgfLp2Rp+Qg+Ffs25cqnUx9l8J9AevELwjhrysIjxa2TH4qXZ4tC5/aqaUktxDQaQuWED+yZMYOjvh8vZb2I8dWyahWVFpEVsitrDk7BLic+JpYdeCr/p9xRDvITdqxR5bpC2Yyr6ipT8Y9iW0GsrxmGvMWXaCg5dSsbMwZtqgVjzRyxs7i+pl91SUxqiqAf8c4AYk1mBbGp3rK2OzdsVoC6YczLB/sJU2j74OAn1xUSmhBxII2RmrBfrWdgx5uh0erSoO9Dn79pG64DsKzpzByM0N1/ffx+7BBzC46ZpAYWkhP138iWXnlpGUl0Q7x3a80fUN+nv21+e5yYTji+HwPMhLA5++MHo+stldBEWmM2fhEY5GpeNkZcLbw9owoYe3qiylKNVQ1b8aJyBMCHGMqle8UiogpaTwUgZZu2IoisvG0M5UWzDVyaVa8+j/Huit7B24+8mpdLi7ioG+sJRzBxII2RVLflYRHq3tuWdyO9xbVhDodTqyd+8m9bvvKAwLx9jDA7ePPsJ2TNnMlfkl+X8F+pT8FDq5dOLjXh/T072nNiSUl64VBD/6vZbnpsVg6Pc60rM7+y6m8O13hwmOzcDF2pT/jPTjkW5emFczF5CiKFUP+B/WZCMak4KIDLJ2xlAUk4Wh7f+WAuFfCfT7EwjZFUN+djFN29jTdUp73FuUP5VRlpaSvWMHqQu+o/DSJUy8vWny+efYjhqJML6RMiGvOI8fL/7IsnPLSCtIo6tbV77s9yVd3bpqO+SmweG5cGyhVhS87Sjo+xo6t0B2hScxd+4hziZk4mFnziej2/NQ56aYVVC/VlGUqqtyxSshhCug/4vlmJSyygXNFa3CVObOaAovZWBgo09q1tWtTgJ9UUEJ5/YncGp3LPnZxXi2tafriGY0uUWgz/p9G6kLFlAUGYlJ8+a4T5+OzbChZVIU5xXnsf7CelaEriC9IJ3uTbozw38GXdz0M8ZyUuDwt1rmyuI8aDcG+r1BqXNbtp1LZO6Gg5y/mo23owVfPeDP6I4emFQzRYSiKP9U1Vk6Y4HpwD5AAN8KId6QUv5Ug21rEIqv5pK5M4aCsDQMLIywHdEMqx5NENU8Y40/H0rQD6uJCztbrUB/dl88p3bHUZBTjJefA11HNsPNt/z0CX8F+vnzKYqKwrRVKzy++RrrIUPKpCjOLc5l3fl1rAhdQUZhBr3ce/FMwDN0dNGXNM5O0gqOHF8CpYXQ/kHo9zolDi3ZcvoK81bvJyIll+bOlnz9cACj/N1VnhtFqQFVHdJ5D+h6/axev9J2N6ACfgVK0vLJ2h1L3qlkhImhVkqwjwcGFSQRq8zViEsc2rCa6FMnsbC1Y8DEKfgPHHp7gX5XHAW5xXi1c6DriEoC/bbtWqCPjMS0ZUs8Zs/GevCgMoE+uyibteFrWRW+iszCTPp69GVqwFQCnAO0HbIS4dA3cHI5lBaD/1jo+xpFds35OSSe+cv3E5OWRxs3a+Y90omh7VVRcEWpSVWNPgZ/G8JJQyuEovxNaWYhWX/Ekns8CWEotOLg/ZpiaFm9ohopMVEc2rCGiBNHMLO2od+EpwgcMhxj08rnnRcXlnJ2XzwhO2MpyC3Gu70jXUb44NbsFoF++3ZS5y+gKCJCC/TffIP1kMFlAn1WURZrwtawKnwV2UXZ9G/an6kBU2/kos+M1+bQB68EWQoB46DPqxTY+JTJXOnf1JaFj3VmUFtXDFSgV5QaV9WAv10IsQNYp//5YeD3mmlS/VSaU0T2vnhyjlwBCZbd3bAZ4IWhTfXmiaclxHH4x7VcOHwQUwtLeo+dQKfh91Yp101JcSmhB65wckcM+VlFeLVzoNtIX1ybVVwYPHv7dlLmz6focgSmLVuUO3STWZjJqrBVrAlfQ05xDnd73s3UgKn4OfppO1yL0YqCh+iTqAY+An1fJd/SkzVHY1h4YC/J2YV09rbnszHtuauVs0popii16JYBXwjRAnCVUr6hL3PYB20M/zCwphbad8fTFZSQfTCBnIMJyOJSLDq5YjPQCyOH6q38zLiayOGN6wg/uA8jU1N63P8wnUeMqVL2ytISHeGHrnBiWwy5GYV4tLan+5T2FV+M1enI3rGDlHnzKLocgUmL5nh8PQvre+4pE+gzCjJYGbaStefXklucy2DvwUzxn0IbhzbaDulRWp6b0+tAGECnx6HPK+RZuLP6SAwLD/xBak4RPX0d+ebhQHo2rzglg6IoNaeyM/xvgHcBpJSbgE0AQogu+vtGVXxowyZLdOQcvkL23jh0eSWYd3DCZrA3xi7Vy+WSlZrMkU0/ELpvNwYGhnQeOZqu9z5QpXz0paU6Lhy5yonfoslOL6BJc1sGPelH04pWxup0ZO/cSeq8eRReuoxJ8+Z4zJqJ9dCh/wj0y0OXs+78OvJL8hniM4Qp/lNoZd9K2yEtQh/o14OBEXSZBL1fJs/cldVHYvh+/17Scovo29KJlwa2pKuPQ7X6RlGUf0dlAd9HSnnm7xullCeEED410qI7nNRJ8k6nkLUjmtKMQkxb2mF7jw8mTauXnTHnWjrHNv/Imd3bAAgYPJxuox/Cyr7y4KjTSS4du8qx36LJSsnHxdua/o+2xtOv/NzvWqDfpQ/0l7TplTNnYDN0KMLwxqyhzMJMVoatZE34GvKK8xjqM5Qp/lNoYd9C2yH1kpai+OwGMDSB7lO1QG/q9I9A/8qglnT2VoFeUe4ElQX8W41LmP+bDakPCi5eI3NbFMWJuRi7W2L/QEvMKliNWulj5eZwfMtGgn/fQmlJMe0HDKbH/Q9j4+RS6bFSJ7kcnMzxX6O4djUPJ08rhj/nj08F2SullGTv2kXq3HkUXryIia8v7jNmaPPobwr02UXZrA5fzarQVWQXZzPEewjPBjx7I9CnRcD+r7RAb2QGPZ+Hni+SZ+rIqsMxLDxwVgV6RbmDVRbwjwshJkspF928UQgxCTh5qwOFEJ5o2TTdAB2wUEo5+39pbF0pSsghc1sUhZczMHQww2GcvspUNWaWFBcVErJtK8d/+YmC3Bza9L6LXmMfxd7NvdJjpZREnUrl2K+RpCXk4uBuydCp7fGtIJumlJLcAwdImT2HgrAwTJo10xZMDR9WJtDnFueyNnwty0OXk1WUxUCvgTwb8CytHfR1ZdOjtOpSp9drZ/Q9X4BeL5FnYq8P9GdUoFeUeqCygP8K8LMQ4lFuBPgugAkwppJjS4DXpJTBQghr4KQQYpeUMux/anEtKkkvIHNHNPmnU7RFUyN9tUVT1Vj9qSst5dy+XRz+aR056Wk0C+xMn/FP4OLjW+mxUkpiQ9M5uiWSlNhs7FwtGDzJjxadK57OmHvkKCmzZ5MfEoJx06Y0+eK/2I4aVSbQX18Zu+zcMjIKM+jftD/PBj57Y9ZNRqwW6E+t1cbouz+jH7q5fkZ/irTcIvq1cublgS3p7F29TzuKotSOWwZ8KWUS0EsIMQDQT7LmNynlH5U9sJQyEX12TSllthAiHPAA7viAX5pbTPYfseQcSUQYCKwHeGJ9V9NqLZqSUnLp6CH+XL+Ka4kJNGnZmuEvvo6nX4cqHZ8YkcmRzRFcuZSBjZMZA59oS6turhWWEsw/dYrk2bPJO3wEI1dX3D78ELsH7i+T6ya/JJ8NFzaw9NxS0gvS6ePRh+cDny87j/7gTAheBUJoF2P7TCPX1JlVR7Qz+nQV6BWl3qlqLp29wN7qPon+Am9H4Gh1H6M2yBIdOUFXyPojFllYimUXN2wGeWFoa1qtx4s5e4qDa1eQFHkJx6Ze3Pf6+zTv0r1KUxLTEnI48ksk0WdSsbAxod+4Vvj1ca+wOHhBWBgps+eQs38/hg4OuL7zNnbjxmFgeqPt19MULz67mNT8VHo26clzgc8R6BKo7ZB1RasuFbxCqx/b6XHo+xq5Zq76QL9XBXpFqcdqPKm4EMIK2Ai8IqXMKuf+KcAUAC8vr5puTrm0SlNpZP4eRUlaAWat7bEd3gxjV8vKDy5HUuRlDq5bQcyZEKwdnbnn2Vfw6zcAA4PK8+dkpuRz7NdILh5LwsTMiB6jffEf4ImxafnHFkZEkDLnW7J37MDAxgbnadNwmPAoBpY32l5UWsSmS5tYdHYRyXnJdHXryoy7ZtDZtbO2Q3aStmDqxFJtZWzHCdD3NfItPFh5OJrvVaBXlAahykXMq/XgQhgDvwI7pJSzKtu/LoqYF13JIfPXSAojMzFyMcduhC9mrat30TEz+SoH167gwuGDmFnb0GPMWAIGD69Svpu8rCJO/B5N6MEEhIEg4O6mdBzijVkFKRmK4uJInTuXzK2/YmBmhsPEJ3CYOLFMzdgSXQlbI7ay4PQCEnMT6eTSiecDn6dbk27aDjkpWq6b40ugtAgCx0O/Nyiw8mTdsVjm7Y0gNadQfzG2lQr0inIHqqki5rfbCAEsAcKrEuxrW2l2EZk7osk7mYSBuZGWrrhbk2pVmirIyeHIzz9wavtWhIEh3cc8TNd778fUovJPCIX5JYTsjOH0njhKSyR+vZvQZXgzrOzLH0YqTkoidd58MjZtQhga4jBxIo6Tn8bI/kYwllKyK2YXc0/NJSozinaO7fiw54c3Co/kpmnZK48thJIC8H8Y+r1BkW0zfjwZx9w/9pGYWUAPXwe+m9CJLmrBlKI0CDU5pNMbeAw4K4Q4pd/2rpSyTnPwyGId2X8mkL03Dlmqw6q3BzZ3e2JgcfvJzUpLijm143eObFxHQV4u7e4aSO+HJ2Dt4FTpsSVFpZzZF0/wjhgKc0to2cWFbqN8sXMtf6VuaWYmaYsXk75yFVKnw37sWByfmYqxy415+1JKgq4EMSdkDmFpYTS3bc43/b/hbq+7tUBfkAVH5kPQXK3wSIcH4a63KLFvzs8hCcz5Yx9x6fl08rJj5kMB9GpR+e+hKEr9UWMBX0r5J1renTtG/oV0MrZEUJpWgJmfozZO73T768eklFw6FsTBNcvJSErEq0Mgd014qkpTLHU6yYUjVzm6JZLcjEK82jnS4z5fnL3KX6mrKyjg2urVpC5ajC4rC5uRI3F+6UVMPD3L7BeSHMLs4NmcTDqJh5UHn/X5jBHNRmjFwYvztZqxB2dBfjq0vRcGvIfOqTVbz1xh9vIDRKbm0t7Dho+fbE9/ldRMURqkRlEJuiS9gIxfIykIS8PI2RynSe2rvUL2ysXz7F+1hCsXw3Fs6sX973yET0CnKgXI2LA0gjZGkJaQg4uPDYOf8qu4QHhJCZmbN5Py7VxKkpKw7NcXl1dfxaxNmzL7XUi/wJyQORyIP4CjmSPvdn+XB1s+iLGhsZaD/sQKbXVsdiI0vxvu/g/SvSM7QpP4es1BLiRl09rVmu8f68wQP1cV6BWlAWvQAV8W68g+EE/W3jiEAJuhPlj38ajWwqmslGT2r1nGxcMHsbC1Y/CUF2jffzAGhpXPvEmNzyFo02XiwtKxcTJjyNPtaNHZpeI0CLt3k/L1NxRFRmIW4I/79K+w7NatzH4xWTHMC5nHtuht2JjY8EqnVxjfZjwWxhag08GZH2HvZ3AtCjy7w/2LkD592HchhZlz/+RcQha+zpbMGd+RkR2aqHz0itIINNiAf/PwjXkHJ2xH+GJkd/vz6YsLCzi+ZSPHf9kIQtDjgXF0vfcBTMwqHwrKuVbI0a2RnD+ciKm5Eb0fbEGHu5piaFz+P5zcY8dImTmL/NOnMfH1xePbOVgPGlTmH8PV3Kt8d/o7Nl/ejImhCZM7TGZi+4nYmNhoc+fP/w5/fArJoeDaAR7ZAC2HEBSZxowFQQTHZuDpYM6MhwIYHahKCSpKY9LgAn5pThEZWyPJP53yPw3fSCm5eORP9q9aSnZaCq179qXfhCerlNysqKCE4B0xnN4dh05KAgd60nmYT4VTLAsuXCB55kxyDxzUVsd+8jF2Y8aUKRCeVZTF4rOLWRO2BolkfJvxTOowCSdz/YXVqAOw52OIPw4OzeHBpeA3hrNXsvlq6TEOXkqlia0Zn4/pwENdmmKsAr2iNDoNJuBLKckLTibzt0h0haXYDPLCur9ntYZvkqMj2btiIfFh53D2bsbwF16jqV/7So/TleoIO5TIsa2R5GcX07KLCz1GN8emggvDxUnJpMyZTeamnzGwtsbl9dewnzABA7MbSUqLSotYd34di84uIqswi1HNR/F84PO4W+mTrSWc1AJ95D6w8YBRcyDwUSLTC5i5/hS/nUnE3sKY90e0ZUIPb8yqWTxdUZT6r0EE/JK0fK79fJnCyxmYeNtg/0DLahUiycvKJGjDas7s3oGplRWDnn6eDgOHVGmFbNz5dP7ccIn0K7m4t7RjxPMtcPUpv6SgLi+PtKXLSFuyBFlSgsPjj+P07DMY2t2oTKWTOrZFbePbkG9JyEmgt3tvpnWediODZVqEFujDNoOFE9zzX+jyFFfzYPYv4Ww4EYepkQEvDWzJ5L7NsDarXk1dRVEajnof8HV5xSTNCQHAbrR+8dRtXoCUOh1n/9jJwbXLKczPI3DoCHo9+GiVygpmpuQTtPEykadSsHEy09IVB5Y/rVGWlpK5+RdSZs+mJDkZ6yFDcHn9NUz+llLiSOIRZp2YRXh6OG0c2vD94O/p5d5LuzMnBfZ/CSeXgaEp3PU29HqBzFIzFuyOYNmhKHRS8lgPb54f0AJn6+rlAVIUpeGp9wHfwMIYu9EtMPW1xagaSc6SoyPZvWQ+iRfP09SvPQOfehYnT+9KjysqKOHkthhO7YnFwNCAHqN9CRjoiVEFQya5R46Q9OVXFIaHYxbgj8c3X2PRqVOZfS6kX+Dr4K85lHAId0t3/tv3vwxvNhwDYQBFuXB4Hhyarc2r7/wE3PU2+aZOLAuK4rt9EWQXljAm0INpg1vh6VC9UouKojRc9T7gA1h2rPxC6t8V5ecR9ONagrdtwczKmmHPv0rbvgMqnYcudZILR69y+OcI8rKKaN3DjZ6jm2NZwQygwogIkqfPIGffPozd3bWSgsOH/2Pmzbch37I1YivWJta83uV1xrUZh6mhKZSWQMhy2PcF5CRB21Ew8AOK7Zuz4UQcs3fvJTm7kIFtXHj9nta0bVL+MJKiKEqDCPi34/oq2b3LF5KTnob/oKH0Gf8E5laV16RNisriwPoLJMdk4+Jjw7BnO+DWrPwi46UZGaR8O5dr69djYG6uXZB97LEy6YrzivNYem4pK0JXoJM6JrabyKQOk7A1tdWmWIb/Cns+gtSL4NkDxq5CenZj27mrTF9+gKjUXLp42zPv0U6qQLiiKJVqVAE/51o6uxfPJ+LEEZy9fBg17W3cW7Wt9LiC3GKObI4g9M8rWNiYMGhiW1p1cyu/rGBpKRk//kjKN7MpzcrCftzDOL3wAkYONwKyTurYGrGVOcFzSM5PZlizYbzS6ZUbM2/ijsHO9yHuKDi1gnFrofVwTsZm8Jl+Ln1rV2uWTuzCgNblL+BSFEX5u0YR8KWUhO7fw76ViygtKqbfo0/SecToSlfJSik5f/gqQZsuU5hbTMAAT7qNaoaJefndlnfiBFc/+5zC8HAsunbF9f33MGvdusw+wUnBfHX8K0LTQung1IGZ/WfeKECSEQu7PoDQTWDlCiO/gY6PEZNRyJdrg/n97FVcrE358oEOPNjZE0O1OlZRlNvQ4AN+VmoyuxbNI/rUSTza+DFk6ss4uHtUelxaQg77110g8XImbr429BsfiLNn+cM+xVevkjx9Blm//YZRkyZ4fD0L66FDy5x5J+Qk8PXJr9kRvQMXCxc+7/M5I3xHaBdkC7O1AiRBc0EYQL83offLXCsx4dvfL7LqSDTGhgZMG9SKyf2aYWHS4F82RVFqQIONHFJKQvftZu+KhUid5O4npxI4ZATC4NYLsYoKSjj+WzSn98Rham7EgMfa0LZn+VM9dYWFpC9bTur330NpKU7PPYvj009jYHFjhkxucS6Lzy5mZehKDIQBzwU8xxPtntDnvCmFkFXwxyfaBdkOY2HQBxRYNGHl4Wi+/eMyuYUlPNzVk2mDWuFiY/aPNiiKolRVgwz4+TnZ7F44l4tHD+Hp14F7nn0ZWxe3So+LC0tn75rzZKcV4NfHnZ6jm2NmVf6CpZw/D3H1448pjo3FevAgXN56C5OmTf+6X0rJr5G/MuvkLFLzUxnpO5KXO72Mm6W+HVEHYcc7cPUsNO0G49aic+/M1jNXmL5jP/HX8hnQ2pl3hrellWvlF5QVRVEq0+ACflzoGX6fN4u8jGv0fWQiXUaNqXSlbEFuMYc2XuZ8UCJ2rhaMea0T7i3tyt23ODmZ5C++IOv3bZj4+OC5ZDFWvXuX2edC+gU+P/o5wcnBdHDqwOwBs/F39tfuTIuAXf8H538FW094YAm0f4CQuAw+WhDEqbgM/JrYsOZpf3qrAiSKovyLGkzAlzodRzb9QNBPa7F3c+eRT2fi6tui0uMigpPZv/4iBTnFdBrqTdcRPuUunpKlpVz74QdSZn2NLCrC6cUXcJw8GYOb6tVmFWUx/9R81p9f8Bf2+QAAE/pJREFUj7WJNR/2/JAxLcfox+lz4OAMbZzeyBTu/g/0fJ6kfMGXG06zKSQBF2tTpj/ozwOdmqp0xYqi/OsaRMAvyMlh27yZRAYfx6/vAAY9/TzGZrce787PLmL/2gtEhKTg5GnFqBcCKqw6lR8aytUPP6Lg7Fkse/XE7f/+DxMfn7/ul1KyNXIrM0/M5FrBNca2HsuLHV+8MZ/+3CZtmmVWAgQ8oo3Tmzmz5M8o5u29TEmp5Ln+zXluQAusTBvES6Ioyh2o3keX/Ows1rz3KtmpqQx86lkChgyvdF569NlU/lh1nsK8YnqM9iVwsBeG5aQL1uXnkzJ7DukrV2Lo4ID7jBnYjCj7+BfSL/DZ0c8ISQ7B39mfBYMW4Ofop92ZHA6/vwHRB8HNHx5chvTsxo7QJD77fT9x6fnc086V94b74eWoUiEoilKz6n3AN7OyplWPPjTv3B2P1rdeRFVcWMqhjZcJPZCAo4cl974UiFPT8hOk5R47RuL7/6E4Nha7cQ/j8uqrGNrcSFuQX5LPd6e/Y0XoCmxMbPi418fc1+I+bfimIEtLcHb0OzCxghEzofOTnE/O5ePFRwmKSKO1qzVrnu6uxukVRak19T7gCyHo98jESvdLis5i19JQMlPyCRzsRfd7m5U7Vq/LzSV55iyurV2LsacnXitWYNm9bHnBoCtBfHL4E+Jz4rm/5f282vnVG8M3ZzZowzc5ydDpcRj4AVmGNsz69TwrD0djY27MJ/e1Y3w3L1VtSlGUWlXvA35lpJSc3RfPoZ8uY2FjwuhXOuLRuvwKWLmHD2tn9VeuYP/4Y7i88kqZOfXpBenMOD6DrZFb8bHxYek9S+nq1lW7My0Cfp0GUfvBozOMX4d078SW01f49LcQUnMKebS7F68PaY2dhUm5z68oilKTGnTAL8ov4Y9V54kITsangyMDJ/qVW2ZQV1BA8vQZXFuzBhMfH7zXrC6Tuvj6Rdnpx6eTU5TDFP8pTPGfomWzLCmCoNmwf7o2+2bETOj8FJdTc/nPoqMcjkwjoKktS57ogn/T8qd6Koqi1IYGG/CvXc3lt/lnyEotoOeY5nQc7FXuatmC8+dJeP11ii5H4PDE4zhPm1amxGByXjIfHf6IA/EHCHAO4IOeH9DSvqV2Z+xR2PoypISD32gY9iV5pk58u/Miiw9GYm5syKej2zO+m5fKe6MoSp1rkAE/Ljyd7QvPYWgkGD0tEPdyiphLnY705StI+fprDO3s8Fy8GKs+NxZQSSn5Leo3/nv0vxSVFvFW17d4pO0j2kXZ/AwtbfGJpdriqfE/QOuh7Ay9ykdbD5CQkc+DnZvy9rA2OFmpilOKotwZGlzADz2YwP51F7F3s2DE8/7YOP6zgHhJWhpX3niT3KAgrAYNpMknn2Bkf+OfQmp+Kp8e+ZQ9sXsIcA7g096f4mPro915Ybt2Vp+bDD2ehwHvklxoxP+tOsn20Ku0cbPmx2d6qvz0iqLccRpUwA/eEcPhnyPwbu/IkEntyk1jnBccQsK0aZRmZOD28UfYPfRQmXn1e2L38FHQR+QU5/Bq51d53O9xDA0MtbP67e/A6bXg0g4eWY9sEsiGE3F89ls4BSU63hzamsl9fTFWs28URbkD1VjAF0IsBUYCyVLK9jX1PNcd/y2KY1ujaNnVlYET2/5jIZWUkmurVpH01XSM3d3xWb8Os7Y35u0XlBQw/fh0NlzcQFuHtvy3739pbtdcu/PiTtj6kjbVst8b0O9NYjKLeUc/p75bMwe+uL8Dvs6VFz1XFEWpKzV5hr8cmAusrMHnALRhnGNbo2jTw40Bj7f9Rx4aXWEhie//h6ytW7EaOBD3/35eZhHVpWuXePPAm1zOuMyT7Z7kxY4vYmxoDAWZsP1dOLUanNvAuLWUuAWy9FAUs3ZdxNjAgM/GtGd8Vy+V+0ZRlDtejQV8KeUBIYRPTT3+dQU5xQRtvIx3e0cGPNbmH4G3JC2N+BdeJD8kBOdXXsFx6pS/hnCklPx48Ue+Ov4VVsZWfD/oe3p59NIOjD0CG5/W8t/0eRX6v03EtWJe/e4wp+MyGNTWlU9Ht8fNVuWoVxSlfqj3Y/hmVsaMfq0Tts7mGPxtGKcwIoK4KVMpSUvDY/ZsbO4Z8td9+SX5fHL4E7ZGbqW3e28+7fMpTuZOUFqiZbXc/6U2A+epneg8urDycDRfbD+PmbEhc8Z3ZJR/E1VLVlGUeqXOA74QYgowBcDLy6taj1Fe6cH80FDiJj0NRkZ4r1qJeYcOf90Xnx3PtH3TuJB+gecCn2Oq/1RtumVGHGyaDLGHtepTI2ZypcCYN5ce48/LqfRv7cyXD/jjqipPKYpSD9V5wJdSLgQWAnTp0kX+G4+ZFxJC3JSpGFhb4b1sGSbe3n/dF3QliDcPvIlOp2PuwLn0a9pPuyN8K/zyvFZ2cMxCCHiYzSEJ/OeXc5TqJJ+P6cD4bp7qrF5RlHqrzgP+v60gLIy4pydj5OSE17KlGLu7/3Xfxosb+eTIJzSzbcbsAbPxsvHShnD++AQOfQPuHeHBpeRaevGfDafYFJxAZ297Zo0NwNvRsg5/K0VRlP9dTU7LXAf0B5yEEPHAB1LKJTX1fABFsbHETp6Cga0NXiuWY+ym1Y+VUvJtyLcsOruI3u69mXHXDKxMrCA3FX56Skt41vlJGPYlYcmFvLDsT6JSc3l5YEteGthSpUVQFKVBqMlZOuNr6rHLo8vLI+6556CkBK9VK/8K9iW6Ej4I+oAtEVt4oOUDvNfjPYwNjCEhGH54DHJT4L55yMBHWX00lk9+DcPO3Jg1T3enV3OVq15RlIajwQzpXP30M4oiIvFashhTX18AinXFvHvwXbZHb+e5wOd4xv8ZbQw+fCtsnAyWzjBpJ3lO7Xlr/Sm2nr7CXa2cmTk2QOXAURSlwWkQAT/32DEyN23CccoULHtp8+iLdcW8uf9Ndsfu5rXOrzGx/UStQEnQXK1AiT5nfVyRFZPnB3EhKZs37mnNs3c1V4uoFEVpkBpEwE9dsAAjNzecnn0G0MbsPwz6kN2xu3mr61tM8JsAOh1sexOOLwK/+2DM9/wZncsL6/5Ep5Msm9iV/q1d6vg3URRFqTn1PuCX5uRSHBOL/biHMTDXMmN+G/ItWyK28Hzg81qwLy3RplyeWQ89X4DBn7DyaCwfbgmlhYsVCx/rgo+TmoWjKErDVu8DvqGVJc137UQWFwPwe+TvLDq7iAdbPchU/6laRapNT0PYLzDgfXR9XuPL7Rf4/kAkg9q6MHtcRyxN6303KIqiVKpBRDphaIgwNCQqM4qPDn9ER5eOvNv9XYTU3Qj293xOQZdneP2HU/x6JpHHe3rzwah2asqloiiNRoMI+AA6qeP9P9/HxNCEr/p9hbEw0lIa64N9XuepTF5xnEOX03hnWBum9PNVq2YVRWlUGkzA33hpI2dSz/B5n89xs3SDPz6F4JXQ93VyOk3lqWXHORGdzqyxAdzfqWldN1dRFKXWNYiAX1xazPenv6ejS0dG+o6E0M1wYDp0fIzsXm8xcekxTsVlMHtcR0YFuFf+gIqiKA1Qg6jFtzt2N0l5SUzxn4JIuQCbn4Om3Si85yumrg7mdFwGc8erYK8oSuPWIM7w98TuwcnciV6uXWHJYDA2R/fQSl7ddJ6giDRmjQ1gWIcmdd1MRVGUOtUgAn54WjidXTtjEDQHEk/D2JV8eSiD384k8s6wNmrMXlEUhQYwpKOTOq7kXsHDxA4OzoK29/J7aTe+PxDJhB5eTOnnW9dNVBRFuSPU+4AvEOx5aA9PpCZDaTGxnd/ijR9P09HLjv8b2U5NvVQURdGr90M6QggcMIIzP6Hr8BCv7MzEyNCA+Y92wsSo3v8/UxRF+dc0jIgYvgWKc/nNdBjBsRl8MMqPJrbmdd0qRVGUO0q9P8MHIGIvOksX3j9mSt+Wdozp6FHXLVIURbnjNIwz/KRzRJi0IbOghLeHtVHj9oqiKOVoEAFfZiVwIsOSwX6utHO3revmKIqi3JHqf8DX6Qhv+xKbC7owoYd3XbdGURTljlX/x/ANDFilG0aYyRV6N3es69YoiqLcser/GT5wJj6Djt72GBk2iF9HURSlRjSICJmYWYCnvZqGqSiKciv1PuDrdJK7WjnTxce+rpuiKIpyR6v3Y/gGBoKvHw6s62YoiqLc8er9Gb6iKIpSNSrgK4qiNBIq4CuKojQSNRrwhRBDhRAXhBCXhRBv1+RzKYqiKLdWYwFfCGEIzAOGAX7AeCGEX009n6IoinJrNXmG3w24LKWMlFIWAeuB+2rw+RRFUZRbqMmA7wHE3fRzvH5bGUKIKUKIE0KIEykpKTXYHEVRlMatJgN+eTmK5T82SLlQStlFStnF2dm5BpujKIrSuNXkwqt4wPOmn5sCV251wMmTJ1OFEDHVfD4nILWax9Yk1a7bo9p1e1S7bk9DbFeV0wQLKf9x0v2vEEIYAReBgUACcBx4REoZWkPPd0JK2aUmHvt/odp1e1S7bo9q1+1p7O2qsTN8KWWJEOIFYAdgCCytqWCvKIqiVK5Gc+lIKX8Hfq/J51AURVGqpiGttF1Y1w2ogGrX7VHtuj2qXbenUberxsbwFUVRlDtLQzrDVxRFUW6h3gf8OyVfjxDCUwixVwgRLoQIFUK8rN/+oRAiQQhxSv81vA7aFi2EOKt//hP6bQ5CiF1CiEv677VaQUYI0fqmPjklhMgSQrxSV/0lhFgqhEgWQpy7aVu5fSQ0c/TvuTNCiE613K7pQojz+uf+WQhhp9/uI4TIv6nvvqvldlX42gkh3tH31wUhxD213K4fbmpTtBDilH57rfTXLWJD7b+/pJT19gtt9k8E4AuYAKcBvzpqSxOgk/62NdqUVD/gQ+D1Ou6naMDpb9u+At7W334b+LKOX8eraPOJ66S/gH5AJ+BcZX0EDAe2oS0u7AEcreV2DQGM9Le/vKldPjfvVwf9Ve5rp/87OA2YAs30f7OGtdWuv90/E/i/2uyvW8SGWn9/1fcz/DsmX4+UMlFKGay/nQ2EU04qiTvIfcAK/e0VwOg6bMtAIEJKWd1Fd/8zKeUBIP1vmyvqo/uAlVJzBLATQjSprXZJKXdKKUv0Px5BW9RYqyror4rcB6yXUhZKKaOAy2h/u7XaLiGEAMYC62riuW/RpopiQ62/v+p7wK9Svp7aJoTwAToCR/WbXtB/NFta20MnehLYKYQ4KYSYot/mKqVMBO0NCbjUQbuuG0fZP8K67q/rKuqjO+l99xTa2eB1zYQQIUKI/UKIvnXQnvJeuzulv/oCSVLKSzdtq9X++ltsqPX3V30P+FXK11ObhBBWwEbgFSllFrAAaA4EAoloHylrW28pZSe0VNXPCyH61UEbyiWEMAHuBX7Ub7oT+qsyd8T7TgjxHlACrNFvSgS8pJQdgVeBtUIIm1psUkWv3R3RX8B4yp5Y1Gp/lRMbKty1nG3/Sn/V94B/2/l6apIQwhjtBV0jpdwEIKVMklKWSil1wCJq6KPsrUgpr+i/JwM/69uQdP1jov57cm23S28YECylTNK3sc776yYV9VGdv++EEE8AI4FHpX7gVz9kkqa/fRJtrLxVbbXpFq/dndBfRsD9wA/Xt9Vmf5UXG6iD91d9D/jHgZZCiGb6M8VxwJa6aIh+fHAJEC6lnHXT9pvH3sYA5/5+bA23y1IIYX39NtoFv3No/fSEfrcngF9qs103KXPWVdf99TcV9dEW4HH9bIoeQOb1j+a1QQgxFHgLuFdKmXfTdmehFR5CCOELtAQia7FdFb12W4BxQghTIUQzfbuO1Va79AYB56WU8dc31FZ/VRQbqIv3V01foa7pL7Qr2hfR/ju/V4ft6IP2sesMcEr/NRxYBZzVb98CNKnldvmizZA4DYRe7yPAEdgDXNJ/d6iDPrMA0gDbm7bVSX+h/dNJBIrRzrAmVdRHaB+55+nfc2eBLrXcrstoY7zX32ff6fd9QP8anwaCgVG13K4KXzvgPX1/XQCG1Wa79NuXA8/8bd9a6a9bxIZaf3+plbaKoiiNRH0f0lEURVGqSAV8RVGURkIFfEVRlEZCBXxFUZRGQgV8RVGURkIFfKVREUKMEUJIIUSbf/ExRwsh/P6tx1OUmqICvtLYjAf+RFuk928ZjZb9UFHuaGoevtJo6HOZXAAGAFuklG30q0N/AGzQajw/CwShrYzsgrZgZqmU8mshRHO0BTHOQB4wGXAAfgUy9V8PACOAZ9Dy3IRJKf/Nfy6KUm01WsRcUe4wo4HtUsqLQoh0fWGJAcAOKeVn+mX2FmjJvzyklO0BhL7ACFrd0WeklJeEEN2B+VLKu4UQW4BfpZQ/6fd/G2gmpSy86VhFqXMq4CuNyXjgG/3t9fqftwJL9cmtNkspTwkhIgFfIcS3wG9oqaWtgF7Aj1pqFEAr6FGeM8AaIcRmYHPN/CqKcvvUkI7SKAghHNFyqySjDdMY6r97o1UkGgG8BEyXUq7UB/h7gIlACvAKcEFK+Y9CFEKI5ZQ9wzdEq7x0L1rOlHbyRsESRakz6qKt0lg8iFZFyFtK6SOl9ASi0AJzspRyEdq4fSchhBNgIKXcCPwHrTxdFhAlhHgI/qo7GqB/7Gy00nUIIQwATynlXuBNwA6wqr1fU1EqpoZ0lMZiPPDF37ZtRMuimCuEKAZygMfRqgst0wdvgHf03x8FFggh3geM0YaFTuu/LxJCvIQ2+2eJEMIWLevh11LKjBr7rRTlNqghHUVRlEZCDekoiqI0EirgK4qiNBIq4CuKojQSKuAriqI0EirgK4qiNBIq4CuKojQSKuAriqI0EirgK4qiNBL/DymPJVeoNidvAAAAAElFTkSuQmCC\n", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHwdJREFUeJzt3Xu4JHV95/H3p6rPzIiAtxkV5sKgTqKgayRHMNE1asQAroDGC6yX1VUnbjQ+Jtn4wKNRw+4+TzTPs14iuZCILK6KqJEdFCWbDe4mu14YFE0YJRnHCyNZB29AZIXTXd/9o6qrq/v0pc6ZU92nZz4vn7arfvWrX31PA9/vr6r6oojAzMwMIJl1AGZmtn64KJiZWclFwczMSi4KZmZWclEwM7OSi4KZmZVcFMzMrOSiYGZmJRcFMzMrtWYdwEpt3rw5du7cOeswzMzmyk033fT9iNgyqd/cFYWdO3eyd+/eWYdhZjZXJH27Tj9fPjIzs5KLgpmZlVwUzMys5KJgZmYlFwUzMyu5KJiZWclFwczMSnP3OQUzs/UmsoBOEFmWP3ei19bptZHl673l5fv12jOir09wv0c/mA3bj2v0b3FRMLOZiKgkwGry6ybRgeQ4mFiX9RmTaMcl6L5xBtpojxonIOsdkyn91H163AYXBTPrN3FWOjC77CbMlc5Ky77t/vHHJegVHTObUiYVkAolSf5cPEgTlKjYVjynCSRCGxKSVPlyq2jrrlf7pdX9k74+3WOWy4PHLPsO36+/T/EsNf5yNVYUJF0O/CvgUEQ8dsh2Ae8GzgHuAV4eEV9qKh47uuWzUsrZ3bKEVklW4xJrLwmPTqxln3HjlMcaFcfgTLh3zGnNSvuSUUtQTXDJiMS6kJBUEmbZZ1TSqyTIvsRaJsQJx6ysj0zQSfOJ9EjS5JnCFcB7gStHbD8b2FU8zgD+uHi2dSSyMUlvSJIbPpvMWDYDrbSPS9DDZ7CrO+ZUiErySyrJbtTsUmghQRvT4TPYymyyPzEOzCYHZ6WjjrnOZqW2/jRWFCLif0naOabLecCVERHA5yU9UNIJEfFPTcU0TctOoSdcH52Y9IYm1gk3sAYTa3vy6X7/DDbLZ9fTkFDOCscmq+qscCEhGZxdDiTQUYm1NztdnlgHjzkssY5M0J6V2pyb5T2FrcBtlfWDRVsjReHe79zFvd/48ehrsEMS5jzceBp/mj0kyVVnpWNnsJUZZatGgh45gx2cLY9I0J6Vmq0LsywKw7LA0FQqaTewG2DHjh2rOth937yLu67/du/IA8lLiaDV0I2nYbPSEddH69zA6hYBEpxMzWxNzbIoHAS2V9a3AbcP6xgRlwGXASwuLq5qDn7sU07k2Cef6BtPZmZjzPITzXuAlyn3JODOJu8nKE1Q90zAzMyGavItqR8GngZslnQQeCuwABARfwJcR/521P3kb0l9RVOxmJlZPU2+++jCCdsDeG1Txzczs5XzF+KZmVnJRcHMzEouCmZmVnJRMDOzkouCmZmVXBTMzKzkomBmZiUXBTMzK7komJlZyUXBzMxKLgpmZlZyUTAzs5KLgpmZlVwUzMys5KJgZmYlFwUzMyu5KJiZWclFwczMSi4KZmZWclEwM7OSi4KZmZVcFMzMrOSiYGZmJRcFMzMruSiYmVnJRcHMzEouCmZmVnJRMDOzkouCmZmVXBTMzKzUaFGQdJakWyXtl3TRkO07JN0g6cuSvirpnCbjMTOz8RorCpJS4FLgbOAU4EJJpwx0ezNwdUQ8AbgA+KOm4jEzs8maPFM4HdgfEQci4j7gKuC8gT4BHF8sPwC4vcF4zMxsgiaLwlbgtsr6waKt6m3ASyQdBK4DfmPYQJJ2S9orae8dd9zRRKxmZkazRUFD2mJg/ULgiojYBpwDfEDSspgi4rKIWIyIxS1btjQQqpmZQbNF4SCwvbK+jeWXh14JXA0QEZ8DNgGbG4zJzMzGaLIo3AjsknSypA3kN5L3DPT5DvDLAJIeQ14UfH3IzGxGGisKEdEGXgdcD3yN/F1Gt0i6RNK5RbffBl4t6SvAh4GXR8TgJSYzM5uSVpODR8R15DeQq21vqSzvA57cZAxmZlafP9FsZmYlFwUzMyu5KJiZWclFwczMSi4KZmZWclEwM7OSi4KZmZVcFMzMrOSiYGZmJRcFMzMruSiYmVnJRcHMzEouCmZmVnJRMDOzkouCmZmVXBTMzKzkomBmZiUXBTMzK7komJlZqVZRkPQ8Sf8o6U5Jd0m6W9JdTQdnZmbT1arZ7x3AcyLia00GY2Zms1X38tH3XBDMzI58dc8U9kr6CHANcG+3MSL+opGozMxsJuoWheOBe4BnVdoCcFEwMzuC1CoKEfGKpgMxM7PZq/vuo22SPiHpkKTvSfq4pG1NB2dmZtNV90bz+4E9wInAVuDaos3MzI4gdYvCloh4f0S0i8cVwJYG4zIzsxmoWxS+L+klktLi8RLgB00GZmZm01e3KPxb4IXA/wX+CXh+0TaWpLMk3Sppv6SLRvR5oaR9km6R9KG6gZuZ2dqr++6j7wDnrmRgSSlwKXAmcBC4UdKeiNhX6bMLuBh4ckT8SNJDV3IMMzNbW2OLgqQ3RsQ7JP0h+ecS+kTE68fsfjqwPyIOFGNdBZwH7Kv0eTVwaUT8qBjv0ArjNzOzNTTpTKH71RZ7VzH2VuC2yvpB4IyBPj8DIOl/Aynwtoj4zCqOZWZma2BsUYiIa4vFeyLio9Vtkl4wYWwNG3LI8XcBTwO2AX8j6bER8eOBY+0GdgPs2LFjwmHNzGy16t5ovrhmW9VBYHtlfRtw+5A+/y0iliLim8Ct5EWiT0RcFhGLEbG4ZYvfCWtm1pRJ9xTOBs4Btkp6T2XT8UB7wtg3ArsknQx8F7gA+NcDfa4BLgSukLSZ/HLSgfrhm5nZWpp0T+F28vsJ5wI3VdrvBn5z3I4R0Zb0OuB68vsFl0fELZIuAfZGxJ5i27Mk7QM6wO9EhD//YGY2I4pY9qai5Z2kBfJ7BI8mvy9wa0Tc13BsQy0uLsbevau5721mdvSSdFNELE7qV/ers88E/hT4BnlxOFnSr0XEpw8jRjMzW2fqFoX/DDw9IvYDSHok8CnARcHM7AhS991Hh7oFoXAA8AfNzMyOMHXPFG6RdB1wNfk9hReQf23F88A/y2lmdqSoWxQ2Ad8DfqlYvwN4MPAc/LOcZmZHDP8cp5mZlWoVheIDaL8B7KzuExEr+uZUMzNb3+peProGeB/5z3BmzYVjZmazVLco/DQi3jO5m5mZzbO6ReHdkt4K/CVwb7cxIr7USFRmZjYTdYvC44CXAs+gd/koinUzMztC1C0KzwUeMavvOzIzs+mo+4nmrwAPbDIQMzObvbpnCg8Dvi7pRvrvKfgtqWZmR5C6ReGtjUZhZmbrQt1PNP/PpgMxM7PZq/uJ5rvJ320EsAFYAH4SEcc3FZiZmU1f3TOF46rrks4HTm8kIjMzm5m67z7qExHX4M8omJlNRZZ1aN93H512u/Fj1b189LzKagIs0rucZGa27kSW0el0yDptsnb+3KksZ50OnXb+3O3T6WtvkxXbx+43Zpxy//ZSub1TGSdrt3sxDq6XMXcg8nT7zFe9lsefeXajr1vddx89p7LcBr4FnLfm0ZjZTEXE6MTWba8kuL5+E/YbnlCHJMnqfmPGGZtg2x0ipvPdnUoS0rRF0kpJ0hZJmpK0WqTpwHqrRZLm7a0Nx+TtxXrSavXWx4zz8Ef9TON/j39Pweww5Yl0BbPGwSTZXW8vDZ0lZpXxlq2PmXV22ku1x+nGGdmUEqmSMomOTYppkUxbKQsbN5EeM9B35L55W1rZlrTSPHlX9hlM5svWW60iMY/YL0lRsqqr8OtW3ctH7wD+I/D/gM8AjwfeEBH/tcHY7AgWEcXp/bjT+YFZ47BEOmq/CafkIy8BjDh9H5eos05nOi+a1Js1jkqoA+utDQskrWN6iWxcQqwm0IkJtVVuHztOuV9/zEdaIj2S1L189KyIeKOk5wIHyX+j+QbARWHKsqwzevY35lS60xkza1zFNc6s3T31rzErHjFLnpaxs8QhSTJttVjYuLE83e87ha+ZjEdeAqjuP3S8XrJdtl+STu01s6NX3aKwUDyfA3w4In4oqaGQmpF1OrSX7lvxDadJN45WfMNpBddJeze7lt9walp1xlgrqbVatDZsnDjbHJ9Qx89Qe0m6/n5KEubt31WzWapbFK6V9HXyy0e/LmkL8NPmwlp7ez/5Cf7mQ1c0eozaN5wqiaveddLWxHHqzV5btS8dOJGaHZ3q3mi+SNLbgbsioiPpJ8zZu4+2n/o4nvriVxzeDSdfJzWzI1zdMwWAxwA7JVX3uXKN42nMCY/6WU541M/OOgwzs3Wt7ruPPgA8ErgZ6L7VIpijomBmZpPVPVNYBE6JmNJdTjMzm4m6F8H/Hnj4SgeXdJakWyXtl3TRmH7PlxSSFld6DDMzWzt1zxQ2A/skfZGav7wmKQUuBc4k/2zDjZL2RMS+gX7HAa8HvrDC2M3MbI3VLQpvW8XYpwP7I+IAgKSryN+xtG+g338A3gH8+1Ucw8zM1lCty0fFL699HTiueHytxq+xbQVuq6wfLNpKkp4AbI+IT44bSNJuSXsl7b3jjjvqhGxmZqtQqyhIeiHwRfKvt3gh8AVJz5+025C28ka1pAR4J/Dbk44fEZdFxGJELG7ZsqVOyGZmtgp1Lx+9CXhiRBwCKD7R/FfAx8bscxDYXlnfBtxeWT8OeCzw2eLTsw8H9kg6NyL21ozLzMzWUN13HyXdglD4QY19bwR2STpZ0gbgAmBPd2NE3BkRmyNiZ0TsBD4PuCCYmc1Q3TOFz0i6Hvhwsf4i4LpxO0REW9LrgOuBFLg8Im6RdAmwNyL2jNvfzMymb2xRkPQo4GER8TvFT3I+hfxeweeAD04aPCKuY6B4RMRbRvR9Ws2YzcysIZMuAb0LuBsgIv4iIn4rIn6TPNG/q+ngzMxsuiYVhZ0R8dXBxuK6/85GIjIzs5mZVBQ2jdl2v7UMxMzMZm/iO4gkvXqwUdIrgZuaCcnMzGZl0ruP3gB8QtKL6RWBRWAD8NwmAzMzs+kbWxQi4nvAL0p6OvkHzQA+FRF/3XhkZmY2dXV/jvMG4IaGYzEzsxnzjwqbmVnJRcHMzEouCmZmVnJRMDOzkouCmZmVXBTMzKzkomBmZiUXBTMzK7komJlZyUXBzMxKLgpmZlZyUTAzs5KLgpmZlVwUzMys5KJgZmYlFwUzMyu5KJiZWclFwczMSi4KZmZWclEwM7OSi4KZmZVcFMzMrNRoUZB0lqRbJe2XdNGQ7b8laZ+kr0r6H5JOajIeMzMbr7GiICkFLgXOBk4BLpR0ykC3LwOLEfEvgI8B72gqHjMzm6zJM4XTgf0RcSAi7gOuAs6rdoiIGyLinmL188C2BuMxM7MJmiwKW4HbKusHi7ZRXgl8usF4zMxsglaDY2tIWwztKL0EWAR+acT23cBugB07dqxVfGZmNqDJM4WDwPbK+jbg9sFOkp4JvAk4NyLuHTZQRFwWEYsRsbhly5ZGgjUzs2aLwo3ALkknS9oAXADsqXaQ9ATgT8kLwqEGYzEzsxoaKwoR0QZeB1wPfA24OiJukXSJpHOLbn8AHAt8VNLNkvaMGM7MzKagyXsKRMR1wHUDbW+pLD+zyeObmdnKNFoUzMzWq4ggy4Ks031kI5aXt3Wq27PBvv3jdIrlqLZlo46zPIZOZf2M5zyCXU98WKOvi4uCmdUSEUQ3mY1Ipp0aSW4wWfZtW2Wy7Ns2MtEPrGdD3wzZCCUiSauPhLSyXLYnvfXWhrRcTlOhVGw6dqHxWF0UzBqUJ6jxs87+meaIZDlq9rrSfSck3c6ERDotEv3JcjB5VpJld721ISFJ077EOm7f3v7Dt00aJ615jCQR0rB36K9PLgq2rnRnop1O1puVrmgGOti3f2a44hnoyCQ6eoZaTazDP5nTALF8BjoxOXZno/19h85gB8dJ1iBZVvp2Z9LVdZsNF4U5t5Lrot1kGcOSZ83romOPMTSBTjjVHxg3pjcZHZ3whp3aF0mrtZCMSXJjZqDLkujh7Dtk3UnU1shRVxQGk2iMm1kOrI9PopMT3uibS6NnneNO7aNIxNOSJBMSU43rohNnoANJuPap/ogZaHfWOTh71Zyd0ptNy1FTFL7837/D5z/xjXV3cylPWEOui45IlmmRcNWXAKdwqu8kanZUOGqKwsN2HsfPPWtHXxLtXcesOXsdkkSXz0J7M11fFzWzeXPUFIUTdz2IE3c9aNZhmJmta/45TjMzK7komJlZyUXBzMxKLgpmZlZyUTAzs5KLgpmZlVwUzMys5KJgZmYlFwUzMyu5KJiZWclFwczMSi4KZmZWclEwM7OSi4KZmZVcFMzMrHTU/J6Cmdk8iAhYWiI6HaLdJtptKJ6T444nPfb+jR7fRcHM5kpkWV+i7EucnQ6x1IZOt71DtJegm2CX2kSnu2/eNrTvUtE2rO9ScZz2Ut+26LSh3NYutne3daC9tLxvu9e/+/eQZSP/9oe/7a086IILGn19XRTMjiAR0UuO1cS5LFmuIJEOJsuBRFmu9yXLztCkOpg4++Kr7jemLzG931kvLSygNEWtFkrTvnVaKWotVLa1UNrK14/ZAK1Wvj1N0UILuttaab6tu77QgrQ7VjH2QN9jTjut8T/VRcGOKtHp5LO2EafnyxNif+Jc21lmZ3hCXDbjHJM8B/4WOp3pv6hJMjFR9iXDbuJstdCmTb2+adpLlEOTaHfcIpGOS6Lltsq+g30r25SmMDjOQtGepkhHz++tuygYUHOGOXFG2R4/uyyXVzib7IxKrDWTZaXvTGaZg8mnmzzL5FhNpr1+yf02jZxVlsmxxgx0JYl0eXJM0cLw+Mr1xO9XOZI0WhQknQW8G0iBP4+I3x/YvhG4Evh54AfAiyLiW03GtFIRkSezIQnx8JNjdxZYSYaTTtHbnfFJdFhCHbWtklBnMsOE4UmmnLkNzCwrM01t3EByzDH9s9LqLLUvIVYT5/jT8/7E2Vred2EgNs8y51ZE0IlO/sg6ZJH1rXeXs6zSXunbjnb+nOXP1X060VnelvXay33qbOu06WRt2lmbc04+h9O3PqnR16WxoiApBS4FzgQOAjdK2hMR+yrdXgn8KCIeJekC4O3Ai5qI585rr+VHH/zQ+BnmsGTZbjcRzmTFTHJcolx2St7dtmFjfxJtpXm/7j6tgW3l+pgZ5tCkOTCjHHWaXuzvhLm2IoIsMjKyMgEFUSayUe2d6OTbiuXutr7EWElSfYkwKxJUZ4mlzlK53M7atDttOt3lrE2WdWhnS7SzDlnR1sk6dLJ25Rjt/mNlbTrFMbPIyKI/ju7f24nibyMjiud8WxR/d+RtRP73FW359vwR6/RfvSSCFEiL5ySgRZAEHH/bIU6/cE6LAnA6sD8iDgBIugo4D6gWhfOAtxXLHwPeK0kRa3+Or1aL5Jj7DUl6C0VCrSbAYduGzzZXlhyr23rJsUz8aaXfOj4l786wus9ZZL3kVElGvUebLO7LlzsZWTsbul+ZoAaTWzFbamdLtDtLxXKbdqeXbNqdJbJOnkDanaVegsmKhJa1yaJDO2vnSaSbiLLiQVYmviwGnsu/qT/mPCl38gRTJqjobesmpeoyQRbFM6OfM1i+rG4beds6TWrDJBEk9BJdGpDQS34J0Oq2BaQsb9tYJMYW1X752C2qyTTfPtiWFvskxThCpMWzQiTky0mI8n8BCQkJgsifVfx/vpSWS91H0n1WgkjzdqUkJKAWSbGspIVI835qkahFkuR9QikogSQFpflzkvLQxzy78X9WTRaFrcBtlfWDwBmj+kREW9KdwEOA7691MH9+17X81S988fAGyYB71yScmQkVCQXIiuWh66O2HWEz+zQC0UtMCfnMrJuwVElQ3fbuckp332K/Ium0KLZ1k1aRhJJivLxvnnBUHLOXgISgTFJ5SsrXe4kIqmsKIal8TiLJ1yvJKimOlKg7Qne5SEqIRHkCS9VtS/IEViathJQ0T3BKSZMUqVX0L7YnCyRKSJJWsZySJC2SJEVpnvCUJChJkYrnvkeeQJXm60llW5LmiTJJ8zGVJCRpq1gWSdoiTVuQiDTpLqekaUpS9E3TPJYkkc9SR2iyKAx7xQfPAOr0QdJuYDfAjh07VhXM8ffbwsPvbPZDH/NAmfoSUTfh9K3TS069WVM1mXXTCoDy5EkvUZUJsJxTdcdT2S8hLRJXN+nkjyTS/D9w5ek0VYtEgiIRJUlKQr6sJC0SUJqPk+QJIi2STaIWadJCSmmlxXKSkiYb8gTWWuglpjRPUknaKtuSIjHlSSopE1SS5rO9bluStFCqIhHmZ3uJRNpqlQkxSVPSJFnXZ4Bm0GxROAhsr6xvA24f0eegpBbwAOCHgwNFxGXAZQCLi4ururT0qnN/j1fxe6vZ1czsqNHktOVGYJekkyVtAC4A9gz02QP8m2L5+cBfN3E/wczM6mnsTKG4R/A64HryS6uXR8Qtki4B9kbEHuB9wAck7Sc/Q2j289tmZjZWo59TiIjrgOsG2t5SWf4p8IImYzAzs/p818vMzEouCmZmVnJRMDOzkouCmZmVXBTMzKykeftYgKQ7gG+vcvfNNPAVGlM0z/HPc+ww3/HPc+ww3/Gvp9hPiogtkzrNXVE4HJL2RsTirONYrXmOf55jh/mOf55jh/mOfx5j9+UjMzMruSiYmVnpaCsKl806gMM0z/HPc+ww3/HPc+ww3/HPXexH1T0FMzMb72g7UzAzszHmuihIOkvSrZL2S7poyPaNkj5SbP+CpJ2VbRcX7bdK+pW6Y67z2L8l6e8k3Sxpb1OxH078kh4i6QZJ/yzpvQP7/HwR/35J71FDP43VUOyfLca8uXg8tInYDzP+MyXdVLzGN0l6RmWf9f7aj4t9Hl770yvxfUXSc+uOOXURMZcP8q/j/gbwCGAD8BXglIE+vw78SbF8AfCRYvmUov9G4ORinLTOmOs19mLbt4DN6/y1vz/wFOA1wHsH9vki8AvkP+T2aeDsOYr9s8DiOn/tnwCcWCw/FvjuHL3242Kfh9f+GKBVLJ8AHKL7y61TyDkreczzmcLpwP6IOBAR9wFXAecN9DkP+C/F8seAXy5mQOcBV0XEvRHxTWB/MV6dMddr7NO06vgj4icR8bfAT6udJZ0AHB8Rn4v8v5wrgfPnIfYpO5z4vxwR3V8/vAXYVMxs5+G1Hxp7AzGOczjx3xMR7aJ9E72fHZ5WzqltnovCVuC2yvrBom1on+IfyJ3AQ8bsW2fMtdBE7JD/i/aXxen17gbiXhbbkBiW9RmIf9yYByeMuRaaiL3r/cXlgd9t6vILaxf/rwJfjoh7mb/Xvhp717p/7SWdIekW4O+A1xTbp5VzapvnojDsH/zgW6lG9Vlp+1prInaAJ0fEacDZwGslPXX1IY51OPEfzphroYnYAV4cEY8D/mXxeOkqYqvjsOOXdCrwduDXVjDmWmgidpiT1z4ivhARpwJPBC6WtKnmmFM1z0XhILC9sr4NuH1UH0kt4AHkP/s5at86Y66FJmKne3odEYeAT9DcZaXDiX/cmNsmjLkWmoidiPhu8Xw38CHW6WsvaRv5vxsvi4hvVPqv+9d+ROxz89pX4v0a8BPyeyPTyjn1zfKGxuE8yG/SHCC/2dq9QXPqQJ/X0n/T5+pi+VT6b9YeIL/hM3HMdRz7/YHjij73B/4PcNZ6e+0r21/O8pu1NwJPonez85x5iL0Yc3OxvEB+Lfk16+21Bx5Y9P/VIeOu69d+VOxz9NqfTO9G80nkiX9znTGn/ZjZgdfoH9I5wD+Q371/U9F2CXBusbwJ+Cj5zdgvAo+o7PumYr9bqbzTYtiY8xA7+bsXvlI8bmky9jWI/1vks6d/Jp8pnVK0LwJ/X4z5XooPV6732MmL8E3AV4vX/t0U7whbT/EDbyafod5ceTx0Hl77UbHP0Wv/0iK+m4EvAeePG3OWD3+i2czMSvN8T8HMzNaYi4KZmZVcFMzMrOSiYGZmJRcFMzMruSiYDZD0XEkh6dFrOOb5kk5Zq/HMmuKiYLbchcDfkn/4aK2cT/6ZBrN1zZ9TMKuQdCz5hwKfDuyJiEcX3yL6EeB48k+g/jvyT4y/j/xDXwFcHhHvlPRI4FJgC3AP8GrgwcAnyb8c7U7yL3R7NvlXcLeBfRGxlgXIbNVasw7AbJ05H/hMRPyDpB9KOo28QFwfEf9JUkr+3fg/B2yNiMcCSHpgsf9l5F+z8I+SzgD+KCKeIWkP8MmI+FjR/yLg5Ii4t7Kv2cy5KJj1uxB4V7F8VbF+LXC5pAXgmoi4WdIB4BGS/hD4FPlXlh8L/CLw0cq3N4/6zv+vAh+UdA1wTTN/itnK+fKRWUHSQ8i/z+gQ+SWhtHg+ifzXsp4NvB74g4i4sigCv0L+BXl3AG8Abo2IE4aMfQX9Zwop8FTgXPLvvjk1ej/CYjYzvtFs1vN84MqIOCkidkbEduCb5Mn7UET8Gfl9hNMkbQaSiPg48LvAaRFxF/BNSS8AUO7xxdh3A8cV7QmwPSJuAN5I/g2gx07vzzQbzZePzHouBH5/oO3jwBXATyQtkX876svIfx3r/UWCB7i4eH4x8MeS3kz+Vc5XkX9z7VXAn0l6Pfm7mt4n6QHkX1X9zoj4cWN/ldkK+PKRmZmVfPnIzMxKLgpmZlZyUTAzs5KLgpmZlVwUzMys5KJgZmYlFwUzMyu5KJiZWen/A8/mMbynHpKmAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] @@ -270,11 +270,52 @@ ], "source": [ "ss = ks_ss()\n", - "plt.plot(ss['a_grid'], ss['c'].T)\n", + "plt.plot(ss['a_grid'][:10], ss['a'][:, :10].T)\n", "plt.xlabel('Assets'), plt.ylabel('Consumption')\n", "plt.show()" ] }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "(7, 500)" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ss['c'].shape" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([1.31353217e-02, 5.84480912e-05, 5.27129890e-05, 3.75346415e-05,\n", + " 3.30050545e-05])" + ] + }, + "execution_count": 26, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ss['D']" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -354,7 +395,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 11, "metadata": {}, "outputs": [], "source": [ @@ -375,7 +416,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 12, "metadata": {}, "outputs": [ { @@ -406,7 +447,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 13, "metadata": {}, "outputs": [ { @@ -441,7 +482,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 14, "metadata": { "scrolled": true }, @@ -474,7 +515,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 15, "metadata": { "scrolled": false }, @@ -529,7 +570,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 17, "metadata": { "scrolled": true }, @@ -553,7 +594,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 18, "metadata": {}, "outputs": [], "source": [ @@ -570,7 +611,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 19, "metadata": {}, "outputs": [], "source": [ @@ -589,7 +630,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 20, "metadata": {}, "outputs": [], "source": [ @@ -606,7 +647,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 21, "metadata": {}, "outputs": [], "source": [ @@ -665,7 +706,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 22, "metadata": {}, "outputs": [], "source": [ @@ -1361,7 +1402,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.8" + "version": "3.6.5" } }, "nbformat": 4, diff --git a/utils.py b/utils.py index 32ccdca..e389b96 100644 --- a/utils.py +++ b/utils.py @@ -8,7 +8,7 @@ '''Part 1: Efficient linear interpolation exploiting monotonicity. Interpolates increasing query points xq against increasing data points x. - + - interpolate_y: (x, xq, y) -> yq get interpolated values of yq at xq @@ -182,7 +182,7 @@ def interpolate_coord_robust_vector(x, xq): ilow = imid else: ihigh = imid - + xqi[iq] = ilow xqpi[iq] = (x[ilow+1] - xq[iq]) / (x[ilow+1] - x[ilow]) @@ -190,7 +190,7 @@ def interpolate_coord_robust_vector(x, xq): '''Part 3: Forward iteration of distribution on grid and related functions. - + - forward_step_1d - forward_step_2d - apply law of motion for distribution to go from D_{t-1} to D_t @@ -354,7 +354,7 @@ def forward_step_transpose_endo_2d(D, x_i, y_i, x_pi, y_pi): alpha = x_pi[iz, ix, iy] beta = y_pi[iz, ix, iy] - Dnew[iz, ix, iy] = (alpha * beta * D[iz, ixp, iyp] + alpha * (1-beta) * D[iz, ixp, iyp+1] + + Dnew[iz, ix, iy] = (alpha * beta * D[iz, ixp, iyp] + alpha * (1-beta) * D[iz, ixp, iyp+1] + (1-alpha) * beta * D[iz, ixp+1, iyp] + (1-alpha) * (1-beta) * D[iz, ixp+1, iyp+1]) return Dnew @@ -363,8 +363,9 @@ def forward_step_transpose_endo_2d(D, x_i, y_i, x_pi, y_pi): '''Part 4: grids and Markov chains''' -def agrid(amax, n, amin=0, pivot=0.25): +def agrid(amax, n, amin=0): """Create grid between amin-pivot and amax+pivot that is equidistant in logs.""" + pivot = np.abs(amin) + 0.25 a_grid = np.geomspace(amin + pivot, amax + pivot, n) - pivot a_grid[0] = amin # make sure *exactly* equal to amin return a_grid @@ -389,11 +390,31 @@ def stationary(Pi, pi_seed=None, tol=1E-11, maxit=10_000): return pi +def mean(x, pi): + """Mean of discretized random variable with support x and probability mass function pi.""" + return np.sum(pi * x) + + def variance(x, pi): """Variance of discretized random variable with support x and probability mass function pi.""" return np.sum(pi * (x - np.sum(pi * x)) ** 2) +def std(x, pi): + """Standard deviation of discretized random variable with support x and probability mass function pi.""" + return np.sqrt(variance(x, pi)) + + +def cov(x, y, pi): + """Covariance of two discretized random variables with supports x and y common probability mass function pi.""" + return np.sum(pi * (x - mean(x, pi)) * (y - mean(y, pi))) + + +def corr(x, y, pi): + """Correlation of two discretized random variables with supports x and y common probability mass function pi.""" + return cov(x, y, pi) / (std(x, pi) * std(y, pi)) + + def markov_tauchen(rho, sigma, N=7, m=3): """Tauchen method discretizing AR(1) s_t = rho*s_(t-1) + eps_t. @@ -452,7 +473,7 @@ def markov_rouwenhorst(rho, sigma, N=7): # invariant distribution and scaling pi = stationary(Pi) s = np.linspace(-1, 1, N) - s *= (sigma / np.sqrt(variance(s, pi))) + s *= (sigma / np.sqrt(variance(s, pi))) y = np.exp(s) / np.sum(pi * np.exp(s)) return y, pi, Pi @@ -539,7 +560,7 @@ def numerical_diff(func, ssinputs_dict, shock_dict, h=1E-4, y_ss_list=None): def numerical_diff_symmetric(func, ssinputs_dict, shock_dict, h=1E-4): """Same as numerical_diff, but differentiate numerically using central (symmetric) difference, i.e. - + f'(xss)*shock = (f(xss + h*shock) - f(xss - h*shock))/(2*h) """ @@ -561,9 +582,9 @@ def numerical_diff_symmetric(func, ssinputs_dict, shock_dict, h=1E-4): def newton_solver(f, x0, y0=None, tol=1E-9, maxcount=100, backtrack_c=0.5, noisy=True): """Simple line search solver for root x satisfying f(x)=0 using Newton direction. - + Backtracks if input invalid or improvement is not at least half the predicted improvement. - + Parameters ---------- f : function, to solve for f(x)=0, input and output are arrays of same length @@ -698,7 +719,7 @@ def printit(it, x, y, **kwargs): def block_sort(block_list, findrequired=False): """Given list of blocks (either blocks themselves or dicts of Jacobians), find a topological sort and also optionally return which outputs must be computed as inputs of later blocks. - + Relies on blocks having 'inputs' and 'outputs' attributes (unless they are dicts of Jacobians, in which case it's inferred) that indicate their aggregate inputs and outputs""" # step 1: map outputs to blocks for topological sort @@ -859,7 +880,7 @@ def __repr__(self): def make_tuple(x): """If not tuple or list, make into tuple with one element. - + Wrapping with this allows user to write, e.g.: "return r" rather than "return (r,)" "policy='a'" rather than "policy=('a',)" @@ -874,11 +895,11 @@ def input_list(f): def output_list(f): """Scans source code of function to detect statement like - + 'return L, Div' - + and reports the list ['L', 'Div']. - + Important to write functions in this way when they will be scanned by output_list, for either SimpleBlock or HetBlock. """ @@ -906,7 +927,7 @@ def extract_dict(savedA, keys, shape): def extract_nested_dict(savedA, keys1, keys2, shape): return {k1: {k2: take_subarray(savedA[k1][k2], shape) for k2 in keys2} for k1 in keys1} - + def take_subarray(A, shape): # verify leading dimensions of A are >= shape