diff --git a/HARK/ConsumptionSaving/ConsIndShockModel.py b/HARK/ConsumptionSaving/ConsIndShockModel.py index 15d4da7f6..659576d4f 100644 --- a/HARK/ConsumptionSaving/ConsIndShockModel.py +++ b/HARK/ConsumptionSaving/ConsIndShockModel.py @@ -366,6 +366,7 @@ def __init__( PermGroFac, BoroCnstArt, MaxKinks, + HyperbolicBeta, ): """ Constructor for a new ConsPerfForesightSolver. @@ -393,6 +394,9 @@ def __init__( additional points will be thrown out. Only relevant in infinite horizon model with artificial borrowing constraint. + HyperbolicBeta: float + Quasi hyperbolic impatience factor in "beta-delta" preferences. + Returns: ---------- None @@ -414,6 +418,7 @@ def __init__( PermGroFac=PermGroFac, BoroCnstArt=BoroCnstArt, MaxKinks=MaxKinks, + HyperbolicBeta=HyperbolicBeta, ) def defUtilityFuncs(self): @@ -612,7 +617,7 @@ def solve(self): The solution to this period's problem. """ self.defUtilityFuncs() - self.DiscFacEff = self.DiscFac * self.LivPrb + self.DiscFacEff = self.DiscFac * self.LivPrb * self.HyperbolicBeta self.makePFcFunc() self.defValueFuncs() solution = ConsumerSolution( @@ -1528,21 +1533,22 @@ def prepareToCalcEndOfPrdvP(self): # Make a dictionary to specify a perfect foresight consumer type init_perfect_foresight = { - 'CRRA': 2.0, # Coefficient of relative risk aversion, - 'Rfree': 1.03, # Interest factor on assets - 'DiscFac': 0.96, # Intertemporal discount factor - 'LivPrb': [0.98], # Survival probability - 'PermGroFac': [1.01], # Permanent income growth factor - 'BoroCnstArt': None, # Artificial borrowing constraint - 'MaxKinks': 400, # Maximum number of grid points to allow in cFunc (should be large) - 'AgentCount': 10000, # Number of agents of this type (only matters for simulation) - 'aNrmInitMean' : 0.0, # Mean of log initial assets (only matters for simulation) - 'aNrmInitStd' : 1.0, # Standard deviation of log initial assets (only for simulation) - 'pLvlInitMean' : 0.0, # Mean of log initial permanent income (only matters for simulation) - 'pLvlInitStd' : 0.0, # Standard deviation of log initial permanent income (only matters for simulation) - 'PermGroFacAgg' : 1.0,# Aggregate permanent income growth factor: portion of PermGroFac attributable to aggregate productivity growth (only matters for simulation) - 'T_age' : None, # Age after which simulated agents are automatically killed - 'T_cycle' : 1 # Number of periods in the cycle for this agent type + "CRRA": 2.0, # Coefficient of relative risk aversion, + "Rfree": 1.03, # Interest factor on assets + "DiscFac": 0.96, # Intertemporal discount factor + "LivPrb": [0.98], # Survival probability + "PermGroFac": [1.01], # Permanent income growth factor + "BoroCnstArt": None, # Artificial borrowing constraint + "MaxKinks": 400, # Maximum number of grid points to allow in cFunc (should be large) + "AgentCount": 10000, # Number of agents of this type (only matters for simulation) + "aNrmInitMean": 0.0, # Mean of log initial assets (only matters for simulation) + "aNrmInitStd": 1.0, # Standard deviation of log initial assets (only for simulation) + "pLvlInitMean": 0.0, # Mean of log initial permanent income (only matters for simulation) + "pLvlInitStd": 0.0, # Standard deviation of log initial permanent income (only matters for simulation) + "PermGroFacAgg": 1.0, # Aggregate permanent income growth factor (only matters for simulation) + "T_age": None, # Age after which simulated agents are automatically killed + "T_cycle": 1, # Number of periods in the cycle for this agent type + "HyperbolicBeta": 1, } @@ -1566,7 +1572,15 @@ class PerfForesightConsumerType(AgentType): MPCmax=1.0, ) time_vary_ = ["LivPrb", "PermGroFac"] - time_inv_ = ["CRRA", "Rfree", "DiscFac", "MaxKinks", "BoroCnstArt"] + time_inv_ = [ + "CRRA", + "Rfree", + "DiscFac", + "MaxKinks", + "BoroCnstArt", + "HyperbolicBeta", + "geometric_solution", + ] poststate_vars_ = ["aNrmNow", "pLvlNow"] shock_vars_ = [] diff --git a/HARK/ConsumptionSaving/ConsIndShockModelDemos/testing_hyperbolic_discounting.ipynb b/HARK/ConsumptionSaving/ConsIndShockModelDemos/testing_hyperbolic_discounting.ipynb new file mode 100644 index 000000000..edf074a5f --- /dev/null +++ b/HARK/ConsumptionSaving/ConsIndShockModelDemos/testing_hyperbolic_discounting.ipynb @@ -0,0 +1,325 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Testing Hyperbolic Discounting" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "I implement hyperbolic discounting by modyfing the `solveOneCycle` function in `core.py`. If a `geometric_solution` attribute is specified, the `solveOneCycle` function uses the solutions specified in this attribute list in its backward induction loop, instead of the current solution. The exponential discount factor is also resized by a factor of the `hyperbolic` discount variable. " + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from HARK.ConsumptionSaving.ConsIndShockModel import PerfForesightConsumerType\n", + "from HARK.utilities import plotFuncsDer, plotFuncs\n", + "\n", + "from copy import copy, deepcopy\n", + "import matplotlib.pyplot as plt\n", + "import warnings; warnings.simplefilter('ignore')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We want to test if this implementation of HARK for hyperbolic agents is accurate. We check the results using a simple 3 period consumption model with Perfect Foresight, for which we have simple algebraic solutions we can use to corroborate our results. The consumption functions have a linear closed form expression, outlined in `quasi_hyperbolic_algebra.pdf`. The slope of the consumption functions are functions of $\\beta$ and $\\delta$ for both the exponential and hyperbolic cases. We solve first for the exponential case. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Exponential Case" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "PFexample = PerfForesightConsumerType()\n", + "PFexample.CRRA = 1.0000000001 # Using a logarithmic utility function\n", + "PFexample.cycles = 1 \n", + "PFexample.T_cycle = 3\n", + "PFexample.T_age = 3\n", + "PFexample.Rfree = 1\n", + "PFexample.PermGroFac = [1, 1, 1]\n", + "PFexample.LivPrb = [1, 1, 1]\n", + "PFexample.aNrmInitStd = 0.0000000000000000000000000000000001" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "unpackcFunc is deprecated and it will soon be removed, please use unpack('cFunc') instead.\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Plot of Consumption Functions\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd1yV5/3/8dfFks1hiDJFRBEVBMUZ98jS7Nlsk9QkNrtm2O+vTdvvN6JxJDHDas1eTZvRpsaY2CaprWlVFJwMFzIUQdn7cM71++M+gqYxKiD3AT7Px6MP8A7nnA+n8Pb2Gp9Laa0RQgjRvbiYXYAQQoiOJ+EuhBDdkIS7EEJ0QxLuQgjRDUm4CyFEN+TWmS8WEhKiY2JiOvMlhRCiy9u2bdtxrXXv83lMp4Z7TEwM6enpnfmSQgjR5SmlDp/vY2RYRgghuiEJdyGE6IYk3IUQohvq1DH3H2K1WiksLKShocHsUs7I09OTyMhI3N3dzS5FCCHOienhXlhYiJ+fHzExMSilzC7nv2itOXHiBIWFhfTv39/scoQQ4pycdVhGKfW6UqpEKbX7lGtBSqkNSql9jo+BbS2goaGB4OBgpwx2AKUUwcHBTv0vCyGE+L5zGXN/E7j0e9eeBv6utR4I/N3x5zZz1mA/ydnrE0KI7ztruGutNwJl37t8FfCW4/O3gKs7uC4hhOjxqpuqWZa+rE2PbetqmT5a66OOz4uBPmf6QqXUXKVUulIqvbS09JxfwNXVleTkZIYNG8YNN9xAXV3dOT/2yJEjXH/99ef89QBTpkyRDVZCCKdg13Y+3fcpsz+dzVt73jr7A35Au5dCauO0jzOe+KG1Xq21TtVap/bufe67Z728vMjMzGT37t14eHjwu9/97pwe19zcTHh4OB999NE5v5YQQjiL3cd3c/u62/nVd78iyi+KD2Z/0KbnaWu4H1NKhQE4Ppa08XnOycSJE9m/fz+1tbXcfffdjB49mpSUFP7yl78A8Oabb3LllVcybdo0pk+fTl5eHsOGDQOMCds5c+aQmJhISkoK33zzDQD19fXcfPPNJCQkcM0111BfX38hvwUhhPhRZQ1lPPPdM9zy+S0U1RTx7IRnefuytxkaPLRNz9fWpZCfAXcCixwf/9LG5zmr5uZmvvjiCy699FKeffZZpk2bxuuvv05FRQWjR49mxowZAGzfvp2dO3cSFBREXl5ey+NfeeUVlFLs2rWL7OxsLr74YnJzc1m5ciXe3t5kZWWxc+dORowYcaG+BSGEOKNmezMf5nzIK5mvUG+t5/Yht3P/8Pvx8/Br1/OeNdyVUh8AU4AQpVQh8AxGqP9RKXUPcBi4sV1V/ID6+nqSk5MB4879nnvuYfz48Xz22WcsXboUMO7K8/PzAZg5cyZBQUH/9Tz/+te/eOihhwAYPHgw/fr1Izc3l40bN/Lwww8DkJSURFJSUkd/C0II8aO2Fm8lbUsa+8r3MTZsLAtGLyDWEtshz33WcNda/+QM/2l6h1RwBifH3L9XCx9//DHx8fGnXd+8eTM+Pj4XshwhhOgwxbXFLEtfxvq89YT7hPP8lOeZHj29Q5ddd6neMpdccgkvvfQSxhwuZGRknPUxEydO5L333gMgNzeX/Px84uPjmTRpEu+//z4Au3fvZufOnReucCGEAJpsTazZtYYr/3wlX+d/zf3D7+fPV/+ZGf1mdPh+GtPbD5yPX/7ylzz66KMkJSVht9vp378/a9eu/dHHzJs3jwceeIDExETc3Nx488036dWrFw888ABz5swhISGBhIQERo4c2UnfhRCiJ9pYuJHFWxaTX53PtKhpPDHqCSL9Ii/Y66mTd8GdITU1VX9/LXlWVhYJCQmdVkNbdZU6hRDOJb8qn8VbF7OxcCMx/jEsGL2A8RHjz+s5lFLbtNap5/OYLnXnLoQQXUWdtY7f7/o9b+15C3cXd34+8ufcmnAr7q6d011Wwl0IITqQ1pov875kafpSjtUd44rYK3hs5GP09j6vI1DbTcJdCCE6SG55Lou2LGJr8VYGBw1myeQlpISmmFKLhLsQQrRTZWMlr2a+yoc5H+Lr4csvx/6S6wZeh6uLq2k1SbgLIUQbnWzw9eL2F6lsquSGQTfwYPKDWDwtZpcm4S6EEG2xs3QnaZvT2H1iNymhKSwYvYCEYOdZUSfhDtx9992sXbuW0NBQdu/effYHCCF6rOP1x3lx+4v8ef+f6e3Vm7SJaczqP8vpDvXpUjtUL5S77rqL9evXm12GEMKJWe1W3tn7Dld8egVrD65lztA5/PWavzI7drbTBTvInTsAkyZNOq2TpBBCnGrz0c0s2rKI/RX7GR8+nqdHP03/gP5ml/WjnCrcf/PXPew9UtWhzzkk3J9nrmhbP2QhRM92tOYoS9KXsOHwBiJ8I3hh6gtMi5rmlHfq3+dU4S6EEM6g0dbIm7vfZM2uNWg085LnMWfoHDzdPM0u7Zw5VbjLHbYQwkxaa74t+Jbntj5HYU0hM/vNZH7qfMJ9w80u7bw5VbgLIYRZ8irzWLR1EZuKNhEbEMvqmasZFz7O7LLaTFbLAD/5yU8YN24cOTk5REZG8tprr5ldkhCik9RZ63h+2/Nc89k1ZJZk8kTqE3x05UddOthB7twB+OCDtp0uLoTourTWrDu0juXpyympL+HKAVfy2MjHCPEKMbu0DiHhLoTocXLKcli4eSHbS7YzJHgIy6YsIzk02eyyOpSEuxCix6hsrOSljJf4U+6f8Pfw55lxz3BN3DWmNvi6UCTchRDdns1u45P9n7Bi+wqqmqq4Kf4mfpb8MwJ6BZhd2gUj4S6E6NYySzJZuHkhWWVZjAgdwS/G/IL4oHizy7rgJNyFEN3S8frjPL/teT478BmhXqEsnriYy/pf1iV2l3YECXchRLditVt5P+t9Vu5YSaOtkXuG3cPcpLl4u3ubXVqnknAHCgoKuOOOOzh27BhKKebOncsjjzxidllCiPP03ZHvWLRlEYcqDzEhYgJPjXqKmIAYs8syhYQ74ObmxrJlyxgxYgTV1dWMHDmSmTNnMmTIELNLE0Kcg6KaIpZuXcrf8v9GpG8kL017icmRk3vMEMwPkXAHwsLCCAsLA8DPz4+EhASKiook3IVwcg3NDbyx+w1e2/0aCsVDKQ9x59A76eXay+zSTOdc4f7F01C8q2Ofs28iXLbonL88Ly+PjIwMxowZ07F1CCE6jNaar/O/Zkn6Eopqirgk5hLmp86nr09fs0tzGs4V7iarqanhuuuu44UXXsDf39/scoQQP+Bg5UEWb1nMd0e+I84Sx2sXv8bosNFml+V0nCvcz+MOu6NZrVauu+46br31Vq699lrT6hBC/LCaphpW7VzFu3vfxcvNi6dGPcVNg2/C3cXd7NKcknOFu0m01txzzz0kJCTw+OOPm12OEOIUWmvWHlzL8m3LOV5/nGviruGREY8Q7BVsdmlOrV3hrpR6DLgX0MAuYI7WuqEjCutMmzZt4p133iExMZHkZKN50MKFC7n88stNrkyIni3rRBYLNy8kszSTYcHDWDF1BYm9E80uq0toc7grpSKAh4EhWut6pdQfgZuBNzuotk4zYcIEtNZmlyGEcKhoqGhp8BXoGchvxv+Gq+OuxkXJERTnqr3DMm6Al1LKCngDR9pfkhCip7LZbXyU+xEvZb5ETVMNtyTcwrzkefh7yAKH89XmcNdaFymllgL5QD3wldb6q+9/nVJqLjAXIDo6uq0vJ4To5rYf207aljSyy7IZ1XcUT49+mkGBg8wuq8tqz7BMIHAV0B+oAP6klLpNa/3uqV+ntV4NrAZITU2VsQ8hxGlK6kpYvm05nx/8nD7efVgyeQmX9LukR+8u7QjtGZaZARzSWpcCKKU+AcYD7/7oo4QQArDarLyb9S6/2/E7rHYrP038Kfcm3tvjGnxdKO0J93xgrFLKG2NYZjqQ3iFVCSG6tU1Fm1i0ZRF5VXlMjpzMk6OeJNpfhm07UnvG3DcrpT4CtgPNQAaO4RchhPghhdWFPLf1Ob4p+IZov2hemf4KkyInmV1Wt9Su1TJa62eAZzqoFtM0NDQwadIkGhsbaW5u5vrrr+c3v/mN2WUJ0W3UN9fz2q7XeGP3G7i6uPLIiEe4Y8gdeLh6mF1atyU7VIFevXrx9ddf4+vri9VqZcKECVx22WWMHTvW7NKE6NK01mw4vIGl6Us5WnuUy/pfxuMjH5cGX51Awh1QSuHr6wsYPWasVqvM1AvRTgcqDpC2JY3NRzczMHAgr094nVF9R5ldVo/hVOG+eMtissuyO/Q5BwcN5qnRT53162w2GyNHjmT//v387Gc/k5a/QrRRdVM1K3es5IOsD/By92LB6AXcGH8jbi5OFTfdnrzbDq6urmRmZlJRUcE111zD7t27GTZsmNllCdFl2LWdvx74K89ve56yhjKuHXgtD494mCDPILNL65GcKtzP5Q77QrNYLEydOpX169dLuAtxjvac2MPCzQvZWbqTpJAkXpn+CkNDhppdVo/mVOFultLSUtzd3bFYLNTX17Nhwwaeesr8v2iEcHZlDWWs2L6CT/Z9QqBnIP970f9y5YArpcGXE5BwB44ePcqdd96JzWbDbrdz4403Mnv2bLPLEsJpNdub+WPOH3k582XqrHXcNuQ2Hhj+AH4efmaXJhwk3IGkpCQyMjLMLkOILiG9OJ20LWnklucyJmwMC0YvYIBlgNllie+RcBdCnJNjtcdYtm0ZXxz6gjCfMJZPWc6M6BmybNhJSbgLIX5Uk62Jt/e+zeqdq7HZbdyXdB/3JN6Dl5uX2aWJH+EU4a61duq//eWUJtFTbSzcyHNbn+Nw1WGmRk3liVFPEOUXZXZZ4hyYHu6enp6cOHGC4OBgpwx4rTUnTpzA09PT7FKE6DQFVQUs3rqYfxT+gxj/GFbOWMmEiAlmlyXOg+nhHhkZSWFhIaWlpWaXckaenp5ERkaaXYYQF1ydtY41u9bw5p43cXdx5/GRj3Nbwm24u7qbXZo4T6aHu7u7O/379ze7DCF6NK01Xx7+kqVbl3Ks7hizYmfx+MjHCfUONbs00Uamh7sQwlz7yvexaMsithRvYXDQYJ6b9Bwj+owwuyzRThLuQvRQVU1VrMxcyQfZH+Dj7sP/G/P/uH7Q9bi6uJpdmugAEu5C9DB2becv+//CC9tfoLyhnOsHXc9DKQ8R6BlodmmiA0m4C9GD7CrdRdqWNHYd30Vy72RWzljJkOAhZpclLgAJdyF6gBP1J3hx+4t8uv9TQrxCWDhhIbNjZzvl8mPRMSTchejGmu3N/CH7D7ya+Sr1zfXcNfQu7ku6D18PX7NLExeYhLsQ3dTW4q0s3LyQ/RX7GRc2jqfHPE1sQKzZZYlOIuEuRDdTXFvM0vSlfJn3JRG+Ebww5QWmRU+TIZgeRsJdiG6i0dbIW3veYs2uNdi1nXnD5zFn2Bw83aR1Rk8k4S5EN/CPgn+waMsiCmsKmRE9g/mj5hPhG2F2WcJEEu5CdGGHqw6zeMti/ln0T/oH9GfVzFWMDx9vdlnCCUi4C9EF1VnrWL1zNW/vfRsPVw/mp87nlsG3SIMv0ULCXYguRGvNF4e+YNm2ZZTUlXDlgCt5dMSj9PbubXZpwslIuAvRReSU5ZC2JY1tx7aREJTAssnLSA5NNrss4aQk3IVwcpWNlbyS+Qof5nyIv4c/vxr3K66Nu1YafIkfJeEuhJOy2W18uv9TVmxfQWVTJTcMuoGHUh4ioFeA2aWJLqBd4a6UsgBrgGGABu7WWv+7IwoToifbUbqDhZsXsvfEXkaEjuAXY35BfFC82WWJLqS9d+4vAuu11tcrpTwA7w6oSYge63j9cZ7f9jyfHfiMUK9QFk1cxOX9L5fdpeK8tTnclVIBwCTgLgCtdRPQ1DFlCdGzWO1WPsj6gJU7VtJga+DuYXczN2kuPu4+Zpcmuqj23Ln3B0qBN5RSw4FtwCNa69pTv0gpNReYCxAdHd2OlxOie/rP0f+waPMiDlQe4KKIi3h61NPEBMSYXZbo4lza8Vg3YASwUmudAtQCT3//i7TWq7XWqVrr1N69ZS2uECcdqTnC498+zk+/+imNtkZWTF3ByukrJdhFh2jPnXshUKi13uz480f8QLgLIU7X0NzAG3ve4PVdrwPwYPKD3DXsLnq59jK5MtGdtDnctdbFSqkCpVS81joHmA7s7bjShOhetNZ8U/ANz219jqKaIi7udzHzU+cT5htmdmnCWTU3Qd4/2/TQ9q6WeQh4z7FS5iAwp53PJ0S3dKjyEIu3LGbTkU0MCBjAmovXMCZsjNllCWfUUAX7N0D257BvAzRWtelp2hXuWutMILU9zyFEd1ZrrWXVjlW8k/UOnq6ePDnqSW4efDPuLtLgS5yi6gjkrIPsdXBoI9it4B0CQ66CwbPgN5ef91PKDlUhLgCtNWsPruX5bc9TWl/K1XFX88iIRwjxCjG7NOEMtIbSbOPuPPtzOLLduB4UC2Pvh8GzIXIUtKPFhIS7EB0suyybhZsXklGSwdDgobww9QWSeieZXZYwm90GBVsge61xl1520LgeMRKm/wriZ0HveOigDWsS7kJ0kIqGCl7OfJk/5f6JAI8Afj3u11wz8BpcVHtWHIsuzVoPB791BPp6qDsOLu4QOxnGPQjxl4P/hZlQl3AXop1sdhsf7/uYFRkrqG6q5ub4m5mXPE8afPVUdWWQu94YbjnwNVjroFcADJxpjJ/HzQBP/wtehoS7EO2QUZJB2uY0ssqySO2TyoIxCxgUOMjsskRnK88zJkNz1sHh70DbwC8ckm8xAr3fBHDz6NSSJNyFaIPSulKWb1vO2oNr6ePdhyWTlnBJzCXS4Kun0BqO7nCscPkcju02rocOgYmPG8Mt4SkdNn7eFhLuQpwHq83Ke1nvsXLHSqx2Kz9N/Cn3Jt6Lt7s0RO32bFY4vMmxwmUdVBWCcoGosXDxszD4cmO1i5OQcBfiHH1X9B1pW9LIq8pjUuQknhr1FNH+0gyvW2ushv1/M8J835fQUAluXjBgGkxdAIMuBR/nXN4q4S7EWRRWF7Jk6xK+LviaKL8oXpn+CpMiJ5ldlrhQqo+1Drcc+gfYmsAryFh7PngWxE4FD+f/l5qEuxBn0NDcwOu7X+f13a/jolx4OOVh7hh6hzT46o5KcyHHsaGoMB3QEBgDo+ca4+dRY8C1a8Vl16pWiE6gtebv+X9nydYlHKk9wqUxl/Lz1J/T16ev2aWJjmK3Q1F66w7RE/uM62HJMPV/jDv00ARTJ0TbS8JdiFMcrDhI2pY0/nP0P8RZ4nj9ktcZ1XeU2WWJjmBtMPq2ZK+FnC+gtgRc3CBmIoy5D+Ivg4BIs6vsMBLuQgA1TTWs3LGS97Pex8vdi6dHP81N8Tfh5iK/Il1afTnkfmUMuez7G1hrwcPX2FAUP8v46GUxu8oLQn5yRY9m13bWHlzL8vTllDWUce3Aa3l4xMMEeQaZXZpoq4oCx4ToWsjbZGwo8u0LSTcak6L9J4Jb9583kXAXPdbeE3tZuHkhO0p3kBiSyMvTX2ZYyDCzyxLnS2tjE9HJ8fPincb1kHi46BFj/Dx8BLj0rB4/Eu6ixylvKGdFxgo+zv2YQM9Afjv+t1wVd5U0+OpKbM2Q/28jzHM+h4p8QBmrWmb+1hhyCYkzu0pTSbiLHqPZ3syfcv/EyxkvU2ut5daEW5mXPA8/Dz+zSxPnoqkW9v/dGHLJXW+Mp7v2ggFTYdITxoYi31Czq3QaEu6iR9h2bBtpm9PIKc9hTN8xPD36aeICe/adXZdQUwq5Xxg7RA9+A80N4GkxgnzwLGOnaC9fs6t0ShLuols7VnuM5duWs+7QOvr69GXZ5GXM7DdTGnw5sxMHWsfPCzYDGgKiYeRdRqBHjwNXOabwbCTcRbfUZGvinb3vsGrnKmx2G3OT5nLPsHukwZczstvhSEbrCUWl2cb1vokw5Wljh2jfxC69ocgMEu6i2/lX0b9YtGURh6sOMyVqCk+OepIovyizyxKnam6EQ/90bPlfBzXFoFwh5iIYOcfosGiRpmztIeEuuo2C6gKe2/oc3xZ8Sz//fqycsZIJERPMLkuc1FAJ+zYYwy37NkBTNbj7QNx0Y/35wJngLfsLOoqEu+jy6pvrWbNrDW/ufhNXF1ceHfEotw+5HQ/Xzj35RvyAyqLWDot5/wK7FXx6w7BrjfHz/pPB3dPsKrslCXfRZWmt+erwVyxNX0pxbTGX97+cx0c+Th+fPmaX1nNpDSVZrR0Wj2QY14PjYNw8Y/15ZCq4uJpbZw8g4S66pP3l+1m0ZRGbizcTHxjPoomLGNlnpNll9Ux2m7Gq5eQKl/JDxvWIVJj+jDHk0lvOlT1fWmsOn6gjo6C8TY+XcBddSnVTNa9mvsoH2R/g4+7D/4z5H64fdL00+OpsTXXGuvPsdcY69LoT4OphDLNc9LCxwsVPWiSfj6oGKzsKKsjIryAjv5zMggrK66xtfj75jRBdgl3b+cv+v/DC9hcobyjnukHX8XDKwwR6BppdWs9Re8LYGZr9ORz4GprroVcADLrYGD+PmwG9ZLfvubDZNbnHqk8L8v2lNWhtrPiM6+3LzCF9SIkOJDnKwpDF5/8aEu7C6e0+vpu0zWnsPL6T4b2H8+qMVxkaPNTssnqGskOO/i3rjF4u2g7+ETDiduPuPGaCbCg6ByXVDWTmV5BRYIT5zsJK6ppsAAR6u5MSHciVw8NJiQ4kKSoAf8/2v6cS7sJplTWU8eL2F/l036cEeQbx7IRnmR07Wxp8XUhaw9FMx/j5OijZY1zvMwwmzjfu0MOGy4aiH9HYbGPPkarT7soLy+sBcHNRDAn354aRkaREB5ISbSE6yPuC7JiWcBdOp9nezIc5H/JKxivUN9dzx5A7uH/4/fh6SA+RC8JmNZYpnrxDryoC5QLR4+GShcYdelB/s6t0SlprCsvr2Z5fboR5QQVZR6postkBCA/wJCU6kLvGx5ASbWFoeACe7p2zUkjCXTiVrcVbSduSxr7yfYwNG8uC0QuItcSaXVb301AF+//m6LD4FTRWgpuXsaFo2v+DgZeAT7DZVTqdmsZmdha0Dq9k5FdworYJAC93VxIjA5gzIYaUKOOuvI+/eWv42x3uSilXIB0o0lrPbn9Joicqri1mWfoy1uetJ9wnnOenPM/06OnS4KsjVRe3big6tBFsTeAdDAlXGMMtsVPAQ3rvnGS3a/aX1rSEeEZ+Bbkl1Wht/PfY3j5MiQ8lJdpCSrSF+D5+uLk6z5BhR9y5PwJkAf4d8Fyih2myNfH23rdZvXM1dm3ngeEPMGfYHLzcvMwurevTGo7nGg25stdBUbpxPbA/jJ5rBHrUGNlQ5HCippHMk0sRC8rZUVBJTWMzAAFe7iRHWbgssa+xgiXSQoC3c08ktyvclVKRwCzgWeDxDqlI9BgbCzeyeMti8qvzmR49nSdGPUGEb4TZZXVtdhsUbm3dUFR2wLgePsIYbhk8G3oP7vETok3NdrKOVhl35Y5Azy+rA8DVRTG4rx9Xp4S3DK/0D/Hpcv+KbO+d+wvAk8AZF7cqpeYCcwGio6XLm4D8qnwWb13MxsKNxPjHsGrGKsZHjDe7rK7L2gAHvzW2/Od8AbWl4OJuHAQ9bp4xIeofbnaVptFac6Sy4ZThlXJ2H6miqdmY9Ozj34uUqEBuHRNNSnQgiREBeHl0/X/NtDnclVKzgRKt9Tal1JQzfZ3WejWwGiA1NVW39fVE11dnreP3u37PW3vewt3FnZ+P/Dm3JtyKu6yTPn91ZbDvK2PIZf/XYK0FDz+js+LgWcZHzwCzqzRFXVMzOwsrW4I8o6CC0upGAHq5uZAUGcCd4/q1LEUMC+ieQ4DtuXO/CLhSKXU54An4K6Xe1Vrf1jGlie5Ca836vPUsTV9KSV0JV8RewWMjH6O3d2+zS+taKvKNsfPstXD4O9A28AuD4Tcb/c9jJoJbL7Or7FR2u+bg8drThldyiquwO24j+4f4MCEuxJj0jApkcJgf7k406XkhtTnctdYLgAUAjjv3+RLs4vtyy3NJ25xG+rF0EoISWDp5KSmhKWaX1TVoDcW7HOvPPzc+B2PMfMKjjg1FKeDSM8IKoLy2iczC1v4rOwoqqGowJj39PN1IjrIwc2ocKdGBDI+yEOTTc9s+yzp3cUFUNlbyauarfJjzIb4evvxy7C+5buB1uMrKjB9nsxp35TnrjLv0ynxAQfRYmPm/RqAHDzC7yk5htdnJKa5uHSsvqODQ8VoAXBTE9/VnVlI4KdEWRkRbiA3xxcWla016XkgdEu5a62+BbzviuUTXZtd2Pt33KS9uf5HKpkpuGHQDDyY/iMXTYnZpzquxBg783bhDz/0SGirAzRNip8LkJ2HQpeDb/Yewik9Oejo2CO0qqqTBakx6hvj2YkS0hRtSI0mJCiQpMgCfXnJv+mPk3REdZmfpThZuXsieE3tICU3hF2N+weCgwWaX5ZxqSoyVLdmfGytdbI3gFWisbBl8OQyYBh4+Zld5wdQ32dh9pPK0DULFVQ0AeLi6MCzCn1tG92vZIBRh8epySxHNJuEu2u14/XFe3P4if97/Z3p79SZtYhqz+s+SX8bvO77fmAzNWQcFWwBtHAI96h7HhqKx4Nr9fiW11uSdqDtleKWcrKPV2ByzntFB3oyJDSI5ykJKdCAJYX70cpPhu/bqfj9JotNY7Vb+kP0HXs18lQZbA3OGzeG+pPvwce++d5znxW6HI9tbd4gezzGuhw2HKQuMQO8ztNttKKqsP+XQiQKjK2KF49AJ315uDI8K4P7JsaREBZIcbSHEt2et8OksEu6iTTYf3cyiLYvYX7Gfi8Iv4qnRT9E/QDoH0txo9G3JdmwoqikG5Wr0PR91L8RfBpYos6vsMM02O7nHasgoaN0gdKDUmPRUCgaF+nHp0L4td+Vxob64yqRnp5BwF+flaM1RlqQvYcPhDUT4RvDi1BeZGjW1Zw/B1FfAvg2ODUV/g6Ya8PA1OiwOnm1sKPLqHidGlVQ1tKwnP3noRL3VOHQi2MeDlGgL146IJDnKQlJkAH4dcOiEaBsJd3FOGm2NvLn7TdbsWoNG87Pkn3HX0AOOPdMAABklSURBVLvwdDOvpampKguNoZacz41e6PZm8AmFxOshfhb0nwTuXfu9abCePHTCWMGSmV9BUYVx6IS7q2JIeAA3jYpq2SAUFSSTns5Ewl38KK013xZ8y3Nbn6OwppCZ/WYyP3U+4b49rFeJ1lCy19GQay0c3WFcDxkE4x407tAjRnbZDUVaa/LL6lq7IuaXs/doFVabMekZYfEiJdrC3RP6kxxlYWi4f6cdOiHaRsJdnFFeZR6Lti5iU9EmYgNi+f3Fv2ds2Fizy+o8tmYo+E/rlv+Kw4CCyFEw4zfGhGjIQLOrbJPqBquj/0rrBqEyx6ET3h6uJEUGcO/EWGOsPMpCqImHToi2kXAX/6XWWsuqnat4Z+87eLp68kTqE/wk4Se4u/SA8dOmOjjwtbFcMecLqC8D114QOxkmPg6DLgO/PmZXeV5sds2+kmrjgGbHCpZ9JTUth07EhfoyfXAoyY7hlUF9fJ3q0AnRNhLuooXWmnWH1rE8fTkl9SVcNeAqHh35KCFeIWaXdmHVHofc9caQy4FvoLne6Kg46FJjU1HcdOh1xq7WTud4TSMZ+RVkOlaw7CiooLbJmPS0eLuTEmVhVqKxbX94lIUArx7wl3YPJOEuAMgpy2Hh5oVsL9nOkOAhLJ+6nOG9h5td1oVz4kBr/5aC/4C2g38kjLjD2CHa7yLoAq2IG5tt7D1SddoJQgVlxqSnm4siIcyf60ZGtixFjAn2lknPHkLCvYerbKzkpYyX+FPunwjwCOCZcc9wTdw13a/Bl9aODUWOM0RLs4zrfRJh0hPG+HnfJKfeUKS1prC8vmXlSkZBOXuKqmiyGf1XwgI8SYm2cPtYo1f5sPDuceiEaBsJ9x7KZrfx8b6PeSnjJaqaqrgp/iZ+lvwzAnp1owMempsg75+tG4qqjxgbivqNh5GLjA1FgTFmV3lGtY3N7CisOGUFSwXHa4xDJzzdXUiKsDDnohiSoywkd+NDJ0TbSLj3QJklmSzcvJCssixG9hnJgtELiA+KN7usjtFQBfs3GIG+bwM0VoG7tzFuHv8rGHQJeAeZXeV/sds1B0prWlauZOSXk3usuuXQidgQHyYNCjFOD4qyEN+35xw6IdpGwr0HOV5/nOe3Pc9nBz4j1DuU5yY9x6Uxl3b9Mdiqo47x88+Nrf92K3iHwJCrjPXnsZPB3bnuastqm8gsKHcMrxjDLNWNxqET/p5uJEcHcsnQviRHW0iOtBDYgw+dEG0j4d4DWO1W3s96n5U7VtJoa+SeYfcwN2ku3u7eZpfWNlpDaU5rh8Wibcb1oFgYe78R6JGjwEnmDZqa7WQXVzlWsBh35Xkn6gDj0InBff25MjmclOhAkqMsxIb4yKETot0k3Lu57458x6ItizhUeYiJERN5avRT9PPvZ3ZZ589uM9rk5nxu3KGXHTSuR4yEab80Ar13vOkTolprjlY2nLYUcVdRJY3NxqRnbz/j0ImbR0eTHGUhMUIOnRAXhvxUdVNFNUUs3bqUv+X/jSi/KF6e9jKToyabXdb5sdYbB1lkr4Wc9VB3HFzcjb4t4x401qD7h5laYl1TM7sKK09bwXKsypj09HBzITEioGX1SnK0hfAAz64/DCa6BAn3bqahuYE3dr/Ba7tfw0W58FDKQ9w59E56uXaRntl1ZadsKPoarHXQyx8GXmysP4+bCZ7+ppRmt2sOnaht6b2SWVBBdnHroRP9gr0ZFxvcMrySEOaPh5tMegpzSLh3E1prvs7/miXpSyiqKeKSmEuYnzqfvj59zS7t7MrzHB0W1xmHQ2sb+IVD8i3G+vN+E8Ct8ycUK+qaWpYhZhYY/6usbz10IjnKwrwpA4ydnpEWguXQCeFEJNy7gYOVB1m8ZTHfHfmOOEscr138GqPDRptd1plpbXRVPLnC5dhu43roEKN/S/zlEJ7SqePnzTY72cXVLcsQMwsqOHjKoRPxffy4PLFvy+lBA3rLoRPCuUm4d2E1TTWs2rmKd/e+i5ebF0+Pfpqb4m/CzcUJ/2+1WeHwptY79MoCUC7GuaEXP2sMuQTFdlo5x6oaTuuIuOuUQydCfD1IjgrkuhGRpERbSIq04CuTnqKLkZ/YLkhrzdqDa1m+bTnH649z7cBreTjlYYK9gs0u7XSN1bD/744NRV9CQyW4ecKAaTDlaaMxl8+Fb0rWYLWxu6iy9UzP/AqOVDYA4OHqwpBwf24eHdWyQSgyUA6dEF2fhHsXs/fEXtI2p5FZmsmw4GGsmLqCxN6JZpfVqvqYo13uOmOli60JvIKMpYrxl8OAqeBx4Q7Q1lpz+ERdy5memQUV7D1SRbNj0jMy0IuRMUHcG2UhJdrCkHB/erk5x3p4ITqShHsXUdFQwYqMFXyU+xGBnoH8dvxvuSruKlyUE6zGKM1tXX9emA5oo2fLqJ8aE6JRY8D1wvyoVTVY2XHK6UGZBRWU1xmTnt4ergyPtDB3UmzLCpbefjLpKXoGCXcnZ7Pb+Cj3I17KfImaphpuTbiVB5IfwN/DnOWAANjtUJTuOHLuczixz7gelgxT/8cYPw8d0uEToja7JvdY9WlBvr/UOHRCKYjr7cvMIX2M4ZVoCwND/WTSU/RYEu5ObPux7aRtSSO7LJtRfUexYPQCBgaadKybtcHo25K91uiwWFsCLm4QMxHG3Gd0WAyI7NCXLKluaOm9kpFfzs7CSuoch04EeruTEh3IlcONbftJUQH4ezp//3UhOouEuxMqqSth+bblfH7wc/p492HJ5CVc0u+Szp/kqy83OitmrzUmRptqwMMXBs6E+FnGRy9Lh7xUY7ONPUeqTrsrLyxvPXRiSLg/N4yMbLkrjw6SQyeE+DES7k7EarPyTtY7rNqxCqvdyk8Tf8q9ifd2boOvioLW9eeHN4G9GXz7QuINxqRo/4ng1r5x65OHTmw/ZSli1pHWQyfCAzxJiQ7krvExpERbGBoegKe7THoKcT4k3J3EpqJNLNqyiLyqPKZETuHJUU8S5R914V9Ya2MTUfY64w69eKdxPSQexj9kBHr4CHBp+8RtTWMzOwtah1cy8is4UdsEgJe7K4mRAcyZEENKlHFX3sffsyO+MyF6tDaHu1IqCngb6ANoYLXW+sWOKqynKKguYMnWJXxT8A39/Pvx6vRXmRg58cK+qK0Z8v/tOKHoc6jIBxREjYaZvzWGXELi2vTUdrtmf2lN6wah/ApyS6rRJw+d6O3DlPhQUqKNpYjxffxwk0MnhOhw7blzbwZ+rrXerpTyA7YppTZorfd2UG3dWn1zPa/teo03dr+Bq4srj4x4hDuG3IGH6wXqodJUa4yb56wzGnPVl4NrL2Pd+cT5xoSob+h5P+2JmsbTDmfeUVBJjePQiQAvd5KjLFyW2NdYihhpIcBbJj2F6AxtDnet9VHgqOPzaqVUFhABSLj/CK01Gw5vYGn6Uo7WHuWy/pfx85E/p49Pn45/sZpSyP3CGHI5+A00N4CnxdgZOniWsVO0l+85P11Ts52so1XGXbkj0PPLjEMnXF0Ug/v6cXVKeMvwSv8QH5n0FMIkHTLmrpSKAVKAzT/w3+YCcwGio6M74uW6rAMVB0jbksbmo5sZFDiIhRMWkto3tWNf5MSB1vXnBZsBDQHRMPIuI9Cjx4Hr2e+etdYcqTyl/0p+ObuPVNHkOHSij38vUqICuXVMNCnRgSRGBODlIZOeQjgLpU8Ohrb1CZTyBf4BPKu1/uTHvjY1NVWnp6e36/W6ouqmalbuWMkHWR/g5e7FQykPccOgGzqmwZfdDkcyWneIlmYb1/smtm7575t41g1FdU3N7CysbAnyjIIKSquNQyd6ubmQFBlAcpSlZSliWIBznUkqRHemlNqmtT6vO8F2pYtSyh34GHjvbMHeE9m1nc8OfMYL216grKHMaPA14mGCPIPa98TNTZC30TEh+gVUHwXlCv3Gw8g5xvh54JmP0rPbNQeP1542vJJTXIWj/Qr9Q3yYEBdiTHpGBTI4zA93mfQUoktpz2oZBbwGZGmtl3dcSd3DnuN7WLhlITtLd5LUO4lXpr/C0JChbX/ChkrHhqLPjY9N1eDuA3HTjeGWgReD9w//pVFe20Rm4en9V6objElPP0/j0ImZU+NIiQ5keJSFIJ/OPxhDCNGx2nPnfhFwO7BLKZXpuPYLrfW69pfVdZU1lLFi+wo+2fcJQZ5B/N9F/8cVA65oW4OvyqLWDUV5/wK7FXx6w7BrHBuKJoP76WvCrTY72UerWw5nziio4NBx49AJFwXxff2ZnRROSrSFEdEWYkN8cZH+K0J0O+1ZLfMvQFLBodnezB9z/sjLmS9Tb63ntiG38cDwB/Dz8Dv3J9EaSrJax8+PZBjXg+Ng3Dxj/XlkKri0Tlweraz/r/4rjY5JzxDfXoyItnBDaiQpUYEkRQbgI4dOCNEjyG96B0gvTidtSxq55bmMCRvDgtELGGAZcG4PttuMVS0nV7iUHzKuR6TC9GeMO/TegwCob7Kx63Bl6115fgXFVa2HTgyL8OfWMf1aNghFWOTQCSF6Kgn3diiuLWZ5+nK+yPuCMJ8wlk9ZzozoGWcP1KY64yCL7M+Ndeh1J8DVwxhmuehhiL8c7duHQ8drjQ1Cm3aTUVBO1tFqbI5Zz+ggb8bEBrWsYEkI85NDJ4QQLSTc26DJ1sTbe99m9c7V2Ow27h9+P3cPuxsvtx9ZHlh7wtgZmrPO2CnaXA+9AmDQxTB4FpXhk8kstRlDLB/lk1mwkwrHoRO+vdwYHhXA/ZNjWw5oDvGVQyeEEGcm4X6eNhZu5Lmtz3G46jDToqbxxKgniPQ7Qx/zskOtE6L5/wZtB/8I7Mm3kh86lX/bBrOtsJaM9eUcKP03YCxHHxTqx6VD+7bclceF+sqhE0KI8yLhfo4KqgpYvHUx/yj8BzH+Mfxuxu+4KOKi079Iazia6Rg/XwclewBoDkkgb/D9/NN1NOuP92Hn5irqrTYgh2AfD1KiLVw7IpLkKAtJkQH4yaETQoh2knA/izprHWt2reHNPW/i7uLO4yMf57aE23A/uYXfZjWWKWZ/btylVxWhlQvFASls6j2PDysT2VoYAIXg7qoYEq65aVRUywahqCCZ9BRCdDwJ9zPQWvNl3pcsTV/KsbpjzI6dzWMjHyPUOxQaqiDrr+icdeicL3FpqsKqerHNfQSfNF/BhuZkyuv9ibB4kRJj4ZeOw5mHhvvLoRNCiE4h4f4D9pXvI21LGluLtzI4aDBLJi8hxSuM+h1/oWzXXwk4+h2u2koF/mxoTuEreyrb3YYzKCyUlOhAFkVZSImyECqHTgghTCLhfoqqpipezXyVP2T/AV93X+6PuYNxRVUEvzEP6vfiBRyz9+Ej+8Xs9Z+AR8xYkvuF8HiUhUF9fOXQCSGE05Bwx2jw9e7uj1i5cwU1zZWMrg/l4WNFDM/5PwD2MICPLXdRH3sp0fEjuSk6kAAvmfQUQjivHhnujc029h6pIiO/gk15/yGnfhWVHicY2tDMMydKGdhURJ7fSDL730PQiKsZ0i+OoTLpKYToQrp9uGutKSyvb+m9kllQQWHREUa7fEd96L/Y6t9AiIuN3x6vY2LwePwvuxq3hEuI8wwwu3QhhGizbhfuNY3N7HS0tz15tufxmkYiKOVyj+382jODncF5rLT4U++iuMszhvuG349v3Exwk12fQojuoUuHu92uOVBa03I4c0Z+BbnHqh2HTmhmBpbwf4E7Ge31b4Kqs9ni2Ytfh/Zlv6uF8UHDeGrC/xEbeI4NvoQQogvpUuFeVtvU0hExs6CCzPwKqhuNQyf8Pd0YEeXPPRGVjLP+m/Dib3CpKoB6RXF0KvNjpvBlzUEifMN5YdQTTIuaJpuHhBDdltOGe1OzneziqtNOD8o7UQcYh04M7uvPlcnhpIZ7ME5n0ufI31G5X0JBBbh5QuxUGic9zltUsibnA+z1FcxLnsecoXPwdJP150KI7s0pwl1rzdHKhtOCfFdR66EToX69SIm2cPPoaFKiLCRaGvE+5Dhy7qtvwdYIXoHG2aGDZ8GAafzjWDqLtiyisKaQGdEzmD9qPhG+EeZ+o0II0UlMCfe6pmZ2FVaS4RhaySgo51hVIwAebi4kRgRw+9h+pEQHkhJtISzAE3XiAGT/Fb5ZBwVbAA2WaBh1D8RfDtHjwNWNw1WHWbzxCf5Z9E9iA2JZPXM148LHmfFtCiGEaTo13Isq6pm14p9kF7ceOtEv2JtxscEtQT64rz8ebi5gt8OR7ZC+yuiweDzHeJKw4TBlgXGH3meo0SMXo8HX6syXeXvv23i4ejA/dT63JNyCu4tsNhJC9DydGu4VdVYCvT2YN2UAKdEWhkdaCD710InmRjj0NWSvhZwvoKYYlCvETIBR9xrDLpao055Ta80Xh75g2bZllNSVcOWAK3ls5GOEeIV05rcmhBBOpVPDfWi4P+/eO+b0i/UVsG+DcSj0vr9BUzV4+ELcdOP80IEzjfH0H5BTlkPaljS2HdtGQlACyyYvIzk0uRO+EyGEcG7mTKhWFhp35tlrjV7o9mbwCYVh1xqB3n8SuJ95RUtlYyWvZL7Chzkf4u/hz6/G/Ypr467F1UXa6QohBHR2uFcXw6rJxmlFAMEDYdyDxvh5RCq4/HhXRZvdxqf7P2XF9hVUNlVy46AbeTDlQQJ6SasAIYQ4VeeHu2sCzPg1xM+C3oPO+aE7SnewcPNC9p7Yy4jQEfxizC+ID4q/YKUKIURX1rnh3nco3LvhvB5yvP44z297ns8OfEaoVyiLJi7i8v6Xy+5SIYT4EZ0b7uexLNFqt/J+1vv8bsfvaLA1cPewu7kv6T683b0vYIFCCNE9OMUO1e/7z9H/kLY5jYOVB5kQMYGnRj1FTECM2WUJIUSX4VThfqTmCEvTl7Lh8AYifSN5adpLTI6cLEMwQghxnpwi3BuaG3hjzxu8vut1AB5KeYg7h95JL1fpry6EEG1harhrrfmm4Bue2/ocRTVFXNzvYuanzifMN8zMsoQQosszLdwPVR5i8ZbFbDqyiThLHGsuXsOYsDFnf6AQQoizale4K6UuBV4EXIE1WutFZ3tMrbWWVTtW8U7WO3i6evLUqKe4afBN0uBLCCE6UJvDXSnlCrwCzAQKga1Kqc+01nvP9JjKxkqu+PQKSutLuSbuGh4Z8QjBXsFtLUEIIcQZtOfOfTSwX2t9EEAp9QfgKuCM4V5YU0iqdyovTH2BpN5J7XhpIYQQP6Y94R4BFJzy50LgvwbNlVJzgbkAof1CeW/We7ioH+8hI4QQon0ueMpqrVdrrVO11qlRIVES7EII0Qnak7RFwKknZ0Q6rgkhhDBZe8J9KzBQKdVfKeUB3Ax81jFlCSGEaI82j7lrrZuVUg8CX2IshXxda72nwyoTQgjRZu1a5661Xges66BahBBCdBCZ3RRCiG5Iwl0IIbohCXchhOiGJNyFEKIbUlrrznsxpaqBnE57QecWAhw3uwgnIe9FK3kvWsl70Spea+13Pg/o7Ja/OVrr1E5+TaeklEqX98Ig70UreS9ayXvRSimVfr6PkWEZIYTohiTchRCiG+rscF/dya/nzOS9aCXvRSt5L1rJe9HqvN+LTp1QFUII0TlkWEYIIbohCXchhOiGOiXclVKXKqVylFL7lVJPd8ZrOiOlVJRS6hul1F6l1B6l1CNm12Q2pZSrUipDKbXW7FrMpJSyKKU+UkplK6WylFLjzK7JLEqpxxy/H7uVUh8opTzNrqmzKKVeV0qVKKV2n3ItSCm1QSm1z/Ex8Fye64KH+ykHaV8GDAF+opQacqFf10k1Az/XWg8BxgI/68HvxUmPAFlmF+EEXgTWa60HA8Ppoe+JUioCeBhI1VoPw2gnfrO5VXWqN4FLv3ftaeDvWuuBwN8dfz6rzrhzbzlIW2vdBJw8SLvH0Vof1Vpvd3xejfELHGFuVeZRSkUCs4A1ZtdiJqVUADAJeA1Aa92kta4wtypTuQFeSik3wBs4YnI9nUZrvREo+97lq4C3HJ+/BVx9Ls/VGeH+Qwdp99hAO0kpFQOkAJvNrcRULwBPAnazCzFZf6AUeMMxRLVGKeVjdlFm0FoXAUuBfOAoUKm1/srcqkzXR2t91PF5MdDnXB4kE6omUEr5Ah8Dj2qtq8yuxwxKqdlAidZ6m9m1OAE3YASwUmudAtRyjv/07m4c48lXYfyFFw74KKVuM7cq56GNtevntH69M8JdDtI+hVLKHSPY39Naf2J2PSa6CLhSKZWHMVQ3TSn1rrklmaYQKNRan/xX3EcYYd8TzQAOaa1LtdZW4BNgvMk1me2YUioMwPGx5Fwe1BnhLgdpOyilFMa4apbWernZ9ZhJa71Aax2ptY7B+Jn4WmvdI+/QtNbFQIFSKt5xaTqw18SSzJQPjFVKeTt+X6bTQyeXT/EZcKfj8zuBv5zLgy54V0g5SPs0FwG3A7uUUpmOa79wnEUreraHgPccN0AHgTkm12MKrfVmpdRHwHaM1WUZ9KA2BEqpD4ApQIhSqhB4BlgE/FEpdQ9wGLjxnJ5L2g8IIUT3IxOqQgjRDUm4CyFENyThLoQQ3ZCEuxBCdEMS7kII0Q1JuAshRDck4S6EEN3Q/wcavMUqXoibhQAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "PFexample.DiscFac = 0.9\n", + "PFexample.solve()\n", + "PFexample.unpackcFunc()\n", + "print(\"Plot of Consumption Functions\")\n", + "for T in range(3):\n", + " x = np.linspace(0,10, 1000,endpoint=True)\n", + " y = PFexample.cFunc[T+1](x)\n", + " plt.plot(x,y, label = T+1)\n", + "plt.xlim([0, 10])\n", + "plt.legend(title='Period')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "According to `quasi_hyperbolic_algebra`:\n", + "\n", + "$$C_3=X_{3}$$\n", + "\n", + "$$C_{2}=\\frac{1}{1+\\delta}X_{2}+\\frac{1}{1+\\delta}$$\n", + "\n", + "$$C_{1}= \\frac{1}{1+\\delta+\\delta^2}X_{1}+\\frac{2}{1+\\delta^2+\\delta^3}$$\n", + "\n", + "Checking if the slopes match the polynomial expressions:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1.0\n", + "0.5263157894710575\n", + "0.36900369003328515\n" + ] + } + ], + "source": [ + "print(PFexample.cFunc[3].derivative(1))\n", + "print(PFexample.cFunc[2].derivative(1))\n", + "print(PFexample.cFunc[1].derivative(1))" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1.0\n", + "0.5263157894736842\n", + "0.36900369003690037\n" + ] + } + ], + "source": [ + "print(1/np.polyval([1],0.9))\n", + "print(1/np.polyval([1,1],0.9))\n", + "print(1/np.polyval([1,1,1],0.9))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "They do." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Hyperbolic Case" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We set a `hyperbolic` factor of 0.7 and assign the previous generated consumption functions as `geometric_solution`." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "unpackcFunc is deprecated and it will soon be removed, please use unpack('cFunc') instead.\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Plot of Consumption Functions\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deVyVdf7//8ebfV8PuIAICCouuCFu4G5mme3LVFNZjU37Zut3mmbm8xvRXLKacnLam6lss5ymdaaaxCZwBxQVVERwgQOyb4dz3r8/LtKmSVNArsPhdb/duolHruu8zkmfXr7O+3q/lNYaIYQQrsXN7AKEEEJ0Pgl3IYRwQRLuQgjhgiTchRDCBUm4CyGEC/LoyiezWCw6Nja2K59SCCG6vc2bN1u11hFnckyXhntsbCybNm3qyqcUQohuTyl14EyPkbaMEEK4IAl3IYRwQRLuQgjhgrq05/5TbDYbJSUlNDU1mV3KSfn4+BAdHY2np6fZpQghxGkxPdxLSkoIDAwkNjYWpZTZ5fwPrTUVFRWUlJQQFxdndjlCCHFafrYto5R6SSlVppTK+8FjYUqpL5RSBW0/hra3gKamJsLDw50y2AGUUoSHhzv1vyyEEOLHTqfn/gpw7o8eexj4l9Y6EfhX28/bzVmD/XvOXp8QQvzYz4a71voboPJHD18IvNr29avARZ1clxBC9Hi1LbUs37S8Xce2d7VML6314bavjwC9TvaNSqkFSqlNSqlN5eXlp/0E7u7ujBw5kmHDhnH55ZfT0NBw2sceOnSIyy677LS/H2Dq1Klyg5UQwik4tIO1BWuZu3Yur+549ecP+AkdXgqpjWkfJ534obVerbVO0VqnRESc/t2zvr6+bNu2jby8PLy8vPjzn/98Wse1trbSt29f3n333dN+LiGEcBZ51jx++fEv+e23v6VfYD/enPtmu87T3nA/qpTqA9D2Y1k7z3Na0tPTKSwspL6+nhtvvJHU1FRGjRrFhx9+CMArr7zCvHnzmD59OjNmzKCoqIhhw4YBxge28+fPZ/jw4YwaNYqvvvoKgMbGRq666iqSkpK4+OKLaWxsPJsvQQghTqmyqZLHv32cq/9xNaV1pfwx7Y+8Nuc1hoYPbdf52rsUch1wPbC47ccP23men9Xa2sonn3zCueeeyx//+EemT5/OSy+9RFVVFampqcycOROALVu2kJOTQ1hYGEVFRcePf/bZZ1FKkZuby65duzjnnHPYs2cPq1atws/Pj/z8fHJychg9evTZeglCCHFSrY5W1uxew7PbnqXR1sgvh/ySX4/4NYFegR0678+Gu1LqTWAqYFFKlQCPY4T620qpm4ADwBUdquInNDY2MnLkSMC4cr/pppuYOHEi69atY9myZYBxVV5cXAzArFmzCAsL+5/zZGZmcueddwIwePBg+vfvz549e/jmm2+46667AEhOTiY5ObmzX4IQQpzSxiMbycjOoOBYAeP7jOeR1EeID4nvlHP/bLhrrX9xkl+a0SkVnMT3Pfcf1cJ7773HoEGD/uvxrKws/P39z2Y5QgjRaY7UH2H5puV8WvQpff378uTUJ5kRM6NTl113q71lZs+ezTPPPIPxGS5s3br1Z49JT0/nb3/7GwB79uyhuLiYQYMGMXnyZN544w0A8vLyyMnJOXuFCyEE0GJv4YXcF5j3wTy+LP6SX4/4NR9c9AEz+8/s9PtpTN9+4Ew89thj3HPPPSQnJ+NwOIiLi+Ojjz465TG33XYbt956K8OHD8fDw4NXXnkFb29vbr31VubPn09SUhJJSUmMGTOmi16FEKIn+qbkG5ZkL6G4tpjp/abzwNgHiA6MPmvPp76/Cu4KKSkp+sdryfPz80lKSuqyGtqru9QphHAuxTXFLNm4hG9KviE2KJZHUh9hYtTEMzqHUmqz1jrlTI7pVlfuQgjRXTTYGvhL7l94dcereLp5cv+Y+7km6Ro83btmd1kJdyGE6ERaaz4r+oxlm5ZxtOEoF8RfwL1j7iXC74xGoHaYhLsQQnSSPcf2sDh7MRuPbGRw2GCWTlnKqMhRptQi4S6EEB1U3VzNc9ueY83uNQR4BfDY+Me4NPFS3N3cTatJwl0IIdrp+w2+ntryFNUt1Vw+8HLuGHkHIT4hZpcm4S6EEO2RU55DRlYGeRV5jIocxSOpj5AU7jwr6iTcgRtvvJGPPvqIyMhI8vLyfv4AIUSPZW208tSWp/ig8AMifCPISM/g/LjznW6oT7e6Q/VsueGGG/j000/NLkMI4cRsDhuv73ydC9ZewEf7PmL+0Pn8/eK/Mzd+rtMFO8iVOwCTJ0/+r50khRDih7IOZ7E4ezGFVYVM7DuRh1MfJi44zuyyTsmpwv33f9/BzkM1nXrOIX2DePyC9u2HLITo2Q7XHWbppqV8ceALogKiWDltJdP7TXfKK/Ufc6pwF0IIZ9Bsb+aVvFd4IfcFNJrbRt7G/KHz8fHwMbu00+ZU4S5X2EIIM2mt+frg1zyx8QlK6kqY1X8WC1MW0jegr9mlnTGnCnchhDBLUXURizcuZkPpBuKD41k9azUT+k4wu6x2k9UywC9+8QsmTJjA7t27iY6O5sUXXzS7JCFEF2mwNfDk5ie5eN3FbCvbxgMpD/DuvHe7dbCDXLkD8Oab7ZsuLoTovrTWfLz/Y1ZsWkFZYxnzBszj3jH3YvG1mF1ap5BwF0L0OLsrd7MoaxFbyrYwJHwIy6cuZ2TkSLPL6lQS7kKIHqO6uZpntj7DO3veIcgriMcnPM7FCRebusHX2SLhLoRweXaHnfcL3+fpLU9T01LDlYOu5PaRtxPsHWx2aWeNhLsQwqVtK9vGoqxF5FfmMzpyNI+Oe5RBYYPMLuusk3AXQrgka6OVJzc/ybq964j0jWRJ+hLmxM3pFneXdgYJdyGES7E5bLyR/wartq+i2d7MTcNuYkHyAvw8/cwurUtJuAMHDx7kuuuu4+jRoyilWLBgAXfffbfZZQkhztC3h75lcfZi9lfvJy0qjYfGPkRscKzZZZlCwh3w8PBg+fLljB49mtraWsaMGcOsWbMYMmSI2aUJIU5DaV0pyzYu45/F/yQ6IJpnpj/DlOgpPaYF81Mk3IE+ffrQp08fAAIDA0lKSqK0tFTCXQgn19TaxMt5L/Ni3osoFHeOupPrh16Pt7u32aWZzrnC/ZOH4Uhu556z93CYs/i0v72oqIitW7cybty4zq1DCNFptNZ8WfwlSzctpbSulNmxs1mYspDe/r3NLs1pOFe4m6yuro5LL72UlStXEhQUZHY5QoifsK96H0uyl/DtoW9JCEngxXNeJLVPqtllOR3nCvczuMLubDabjUsvvZRrrrmGSy65xLQ6hBA/ra6ljudznuevO/+Kr4cvD419iCsHX4mnm6fZpTkl5wp3k2ituemmm0hKSuK+++4zuxwhxA9orflo30es2LwCa6OVixMu5u7RdxPuG252aU6tQ+GulLoXuBnQQC4wX2vd1BmFdaUNGzbw+uuvM3z4cEaONDYPWrRoEeedd57JlQnRs+VX5LMoaxHbyrcxLHwYT097muERw80uq1tod7grpaKAu4AhWutGpdTbwFXAK51UW5dJS0tDa212GUKINlVNVcc3+Ar1CeX3E3/PRQkX4aZkBMXp6mhbxgPwVUrZAD/gUMdLEkL0VHaHnXf3vMsz256hrqWOq5Ou5raRtxHkJQsczlS7w11rXaqUWgYUA43A51rrz3/8fUqpBcACgJiYmPY+nRDCxW05uoWM7Ax2Ve5ibO+xPJz6MANDB5pdVrfVkbZMKHAhEAdUAe8opa7VWv/1h9+ntV4NrAZISUmR3ocQ4r+UNZSxYvMK/rHvH/Ty68XSKUuZ3X92j767tDN0pC0zE9ivtS4HUEq9D0wE/nrKo4QQArDZbfw1/6/8efufsTls/Gr4r7h5+M09boOvs6Uj4V4MjFdK+WG0ZWYAmzqlKiGES9tQuoHF2YspqiliSvQUHhz7IDFB0rbtTB3puWcppd4FtgCtwFba2i9CCPFTSmpLeGLjE3x18CtiAmN4dsazTI6ebHZZLqlDq2W01o8Dj3dSLaZpampi8uTJNDc309raymWXXcbvf/97s8sSwmU0tjbyYu6LvJz3Mu5u7tw9+m6uG3IdXu5eZpfmsuQOVcDb25svv/ySgIAAbDYbaWlpzJkzh/Hjx5tdmhDdmtaaLw58wbJNyzhcf5g5cXO4b8x9ssFXF5BwB5RSBAQEAMYeMzabTT6pF6KD9lbtJSM7g6zDWSSGJvJS2kuM7T3W7LJ6DKcK9yXZS9hVuatTzzk4bDAPpT70s99nt9sZM2YMhYWF3H777bLlrxDtVNtSy6rtq3gz/018PX15JPURrhh0BR5uThU3Lk/e7Tbu7u5s27aNqqoqLr74YvLy8hg2bJjZZQnRbTi0g7/v/TtPbn6SyqZKLkm8hLtG30WYT5jZpfVIThXup3OFfbaFhIQwbdo0Pv30Uwl3IU7TjoodLMpaRE55DsmWZJ6d8SxDLUPNLqtHc6pwN0t5eTmenp6EhITQ2NjIF198wUMPmf8XjRDOrrKpkqe3PM37Be8T6hPK/036P+YNmCcbfDkBCXfg8OHDXH/99djtdhwOB1dccQVz5841uywhnFaro5W3d7/Nn7b9iQZbA9cOuZZbR9xKoFeg2aWJNhLuQHJyMlu3bjW7DCG6hU1HNpGRncGeY3sY12ccj6Q+woCQAWaXJX5Ewl0IcVqO1h9l+eblfLL/E/r492HF1BXMjJkpy4adlIS7EOKUWuwtvLbzNVbnrMbusHNL8i3cNPwmfD18zS5NnIJThLvW2qn/9pcpTaKn+qbkG57Y+AQHag4wrd80Hhj7AP0C+5ldljgNpoe7j48PFRUVhIeHO2XAa62pqKjAx8fH7FKE6DIHaw6yZOMS/l3yb2KDYlk1cxVpUWlmlyXOgOnhHh0dTUlJCeXl5WaXclI+Pj5ER0ebXYYQZ12DrYEXcl/glR2v4OnmyX1j7uPapGvxdPc0uzRxhkwPd09PT+Li4swuQ4geTWvNZwc+Y9nGZRxtOMr58edz35j7iPSLNLs00U6mh7sQwlwFxwpYnL2Y7CPZDA4bzBOTn2B0r9FmlyU6SMJdiB6qpqWGVdtW8eauN/H39Oc3437DZQMvw93N3ezSRCeQcBeih3FoBx8WfsjKLSs51nSMywZexp2j7iTUJ9Ts0kQnknAXogfJLc8lIzuDXGsuIyNGsmrmKoaEDzG7LHEyDgcczGrXoRLuQvQAFY0VPLXlKdYWrsXia2FR2iLmxs91yuXHArAWQM4a47+q4nadQsJdCBfW6mjlrV1v8dy252hsbeSGoTdwS/ItBHgFmF2a+LF6K+S9B9vfgkNbQLlB/FSY9hv4/VVnfDoJdyFc1MYjG1mUtYjCqkIm9JnAw+MeJj443uyyxA/ZGmH3x7B9DRT+E7Qdeg+Hc/4Iwy+DwO9nzUq4C9HjHak/wrJNy/is6DOiAqJYOXUl02OmSwvGWTgccCDTaLnsXAfNNRAUBRPvhOQroVfnfAYi4S6Ei2i2N/Pqjld5IfcFHNrBbSNuY/6w+fh4yNYZTqFsF+S8BTnvQE0JeAXCkAsh+QqITYNOXoIq4S6EC/j3wX+zOHsxJXUlzIyZycKxC4kKiDK7LFF7FPLeNfroR3JAuUPCDJj1exh0Hnj5nbWnlnAXohs7UHOAJdlLWF+6nrjgOJ6f9TwT+040u6yeraUedv3DCPR9X4F2QN9RcO4SGHYJBHTNlg4S7kJ0Qw22BlbnrOa1na/h5e7FwpSFXD34atngyywOO+z/N+S8Dfl/h5Y6CI6BtPuMPnrEwC4vScJdiG5Ea80n+z9h+ebllDWUMW/APO4ZfQ8RfhFml9YzHckz+ui570LtYfAOhmGXGoEeMwHczBsULuEuRDexu3I3GdkZbD66maSwJJZPWc7IyJFml9Xz1ByC3HeM5YtlO8DNAxLPgeTFMPBc8HSOD7Al3IVwctXN1Ty77VnW7F5DkFcQv53wWy5JuEQ2+OpKzbVGu2X7W7D/G0BD9Fg4bxkMvQT8w82u8H9IuAvhpOwOO2sL1/L0lqepbqnm8oGXc+eoOwn2Dja7tJ7B3mp8IJqzBvI/gtZGCI2FKQ8ZyxfDB5hd4Sl1KNyVUiHAC8AwQAM3aq3/0xmFCdGTbS/fzqKsReys2MnoyNE8Ou5RBoUNMrss16c1HN5uBHruu1BfBr6hMPJqo4/eLxW6yc1gHb1yfwr4VGt9mVLKCzh7izaF6AGsjVae3Pwk6/auI9I3ksXpizkv7jy5u/RsqzoIuW8bfXTrbnD3goGzIfkqo5/u4WV2hWes3eGulAoGJgM3AGitW4CWzilLiJ7F5rDxZv6brNq+iiZ7EzcOu5EFyQvw9/Q3uzTX1VQNOz80Av1ApvFYzASYuxKGXmRcsXdjHblyjwPKgZeVUiOAzcDdWuv6H36TUmoBsAAgJiamA08nhGv67vB3LM5azN7qvUyKmsTDYx8mNjjW7LJck91mbNCVswZ2fwKtTRCeYOy8mHy50VN3EUpr3b4DlUoBvgMmaa2zlFJPATVa68dOdkxKSoretGlT+yoVwsUcqjvEsk3L+OLAF0QHRPPg2AeZ2m+qtGA6m9ZQusVYj573HjRUgF84DLvM6KNHjXb6PrpSarPWOuVMjunIlXsJUKK1/n5MyLvAwx04nxA9QlNrEy/veJmXcl8C4I6Rd3DDsBvwdvc2uTIXc6zIuGM0Zw1UFIK7Nww+z+ijJ8wAF7+bt93hrrU+opQ6qJQapLXeDcwAdnZeaUK4Fq01Xx38iic2PkFpXSnn9D+HhSkL6RPQx+zSXEfjMdix1uijH/zOeCw2HSbdA0PmgU/PWUba0dUydwJ/a1spsw+Y3/GShHA9+6v3syR7CRsObWBA8ABeOOcFxvUZZ3ZZrqG1GQo+N67Q93wG9haIGAwzHofhl0NIP7MrNEWHwl1rvQ04oz6QED1Jva2e57c/z+v5r+Pj7sODYx/kqsFX4enm2i2Bs05rOJjd1kd/H5qqwD8Sxt5s9NH7jHD6PvrZJneoCnEWaK35aN9HPLn5Scoby7ko4SLuHn03Fl+L2aV1bxV7TwyOPlYEHr6QNNfoo8dPBXeJtO/JOyFEJ9tVuYtFWYvYWraVoeFDWTltJckRyWaX1X3VV8CO9419XUo3AQrip8CUh41g9w40u0KnJOEuRCepaqriT9v+xDt73iHYK5jfTfgdFydejJsyb9vXbsvWBHs+Na7QCz4HRyv0Ggaz/s8YHB3U1+wKnZ6EuxAdZHfYea/gPZ7e+jS1LbVcNegqbht5m2zwdaYcDij+j9FH3/EhNFdDYB8Yf6vRduk9zOwKuxUJdyE6YGvZVjKyMsivzCelVwqPjHuEgaFdP3WnWyvfc2JwdHUxePobyxaTr4S4yZ0+OLo7aLLZ2XzgGOsLrGwotLbrHBLuQrRDeUM5Kzav4KN9H9HLrxdLJy9lduxsubv0dNWVnxgcfXgbKDcYMB1m/Na40cirZ+2p43Bo8o/UkFlgJbPQSvb+SppbHXi4KUbHtG+PGwl3Ic6AzW7jb/l/Y9X2VdgcNn41/FfcPPxm/DxlQ9Sf1dIAuz82+uiF/wJtN5Yszs4wRtMF9jK7wi51qKqRzEIrmW1X5xX1xr6LiZEB/CI1hvREC+Piwwnw9kDdeubnl3AX4jR9W/otGdkZFNUUMTl6Mg+NfYiYINkM75QcdijKNAJ95zpoqYWgaJh0l9FHjxxsdoVdprbJxnf7KsksKGd9oZV95cYei5YAb9ITLaQlRpCWYKF3cOeM6ZNwF+JnlNSWsHTjUr48+CX9Avvx7IxnmRw92eyynNvRnSf66LWHwDsIhl5oBHr/SaYOju4qNruD7QerWN/Watl2sAq7Q+Pj6ca4uHCuTo0hLdHCoF6BZ6WdJ+EuxEk0tTbxUt5LvJT3Em7KjbtG3cV1Q6+TDb5OpvbIicHRR3ONwdEJM2H2H2HQHPD0NbvCs0przd7yejILyskstPLdvkrqmltRCpKjgvn1lHjSEiIY3T8Eb4+z/yGxhLsQP6K15l/F/2LpxqUcqj/EubHncn/K/fT27212ac6nuQ52/cO4St/3NWgHRI2BOUth2CXg79p35FrrmtnQ1jfPLLRyuLoJgJgwP+aN7EtagoWJA8IJ8ev6SU4S7kL8wL6qfWRkZ/Dd4e9ICEngpdkvMbb3WLPLci4OuxHk3w+OttVDSAyk328sX7Qkml3hWdNks5O9v5LMQivrC6zkH64BINjXk4kDwrljuoX0hAhiws3/gF3CXQigrqWOVdtX8Ub+G/h6+vJw6sNcOehKPNzkjwhgbNR1JLdtcPQ7UHfU2D43+XKjjx4z3iU36nI4NDsO1bC+sJzMAiubDhyjpdWBp7tiTP9QHpg9iLQEC8OignF3c67XL79zRY/m0A4+2vcRKzatoLKpkksSL+Gu0XcR5hNmdmnOobr0xODo8nxw82wbHH2lMTjas3NWdjiTg5UNx5cofrvXyrEGGwCDewdy3fj+pCVaSI0Lw8/LuePTuasT4izaWbGTRVmL2F6+neGW4fxpxp8YZpFb3Gmqgfy/G330/esBDf3GwfkrYOjF4Odaf/FVN9r4z17r8UAvqmgAoFeQN9MH9yItMZxJCRYiA7vXX2QS7qLHOdZ0jKe3Ps17e94j1CeUP0z8AxcmXNizN/iy22DvV0ag7/oYWhshLB6mPgzJVxhfu4iWVgdbi48d75vnlFTh0ODn5c74+HCumxBLeqKFhMiAbn3HsYS76DFaHa28s+cd/rT1T9Tb6rkm6RpuG3kbgV49dMtYreHQ1rY++rvQYAXfUBh1jdFHj05xiT661pqCsjpjvXlBOVn7K2loseOmYES/EO6YlkBaYgQj+4Xg5eE6f8FLuIseYfPRzWRkZbD72G7G9R7Hw6kPkxCaYHZZ5qgqbht48TZY9xiDowed2zY4eiZ4dP2yvc5WVtN0vM2SWWilrLYZgDiLP5eOjiYt0cL4+HCCfV13IpaEu3BpR+uPsmLzCj7e/zG9/XuzfMpyZvWf1a3/ud0ujVWw80Mj1A9sMB7rPwkm3A5DLgLfEHPr66CGllay9lUeD/TdR2sBCPXzZFKChbQEC2mJFqJDzV+i2FUk3IVLarG38PrO13k+53nsDjsLkhdw07CbetYGX60tUPhPo4+++1OwN0N4Ikz/DQy/AkL7m11hu9kdmtzSamOflgIrW4qPYbNrvDzcGBsbykWjBpOeaGFInyDcnGyJYleRcBcuJ7M0k8XZizlQc4Cp/aby4NgH6RfYz+yyuobWULLpxODoxkrws0DKfGP5Yt9R3baPfqCivq1vbixRrGlqBWBInyBunBRHWqKFsbFh+Hj2vP3ff4qEu3AZB2sP8sTGJ/j64Nf0D+rPqpmrSItKM7usrlG5z+ih56wxvvbwgcHnG330AdPAvfv1lqsaWthQWEFmobFXy8HKRgD6Bvtw7rDepCVGMHFAOJYA2evnp0i4i26vsbWRF3Jf4JW8V3B3c+ee0ffwyyG/xMu9+38weEoNlbBjrRHoB7MABXHpxjYASfPAJ8jsCs9Ic6udzUXGEsXMQiu5pdVoDYHeHowfEM7NafGkJVqIt/j3vM9M2kHCXXRbWms+P/A5yzYt40j9Ec6LO4/7xtxHL38XHvrQ2gx7PjMCfc9n4LBBRBLM/B0MvxyCo82u8LRprck/XMuGQivrC61k76+gyWZMHxoVE8LdMxJJT7QwIjoED3fXWaLYVSTcRbdUeKyQxdmLyTqSxaDQQSxOX8yYXmPMLuvs0BqKv2sbHL0WmqohoBeMu8Xoo/ce3m366Eeqm1jftiXuhkIr1jpj+tCACH+uGhtDWoKFcfFhBPp0vzaSs5FwF91KbUstz217jjd3vYm/pz//b9z/47KBl7nmBl/WwraBF2uMtemefpB0Qdvg6Cng7vyvua65le/2VhxvtRSW1QFgCfD6ryWKfYJde693Mzj/7w4hMDb4+rDwQ1ZuWcmxpmNcOvBS7hp1F6E+7Rse7LTqrZD3nhHopZuNwdHxU2Ha/4PBc8E7wOwKT6nV7mB7SRWZBcYHoVuLq2htmz6UGhfOFSnRpCVEMLh3YI9dothVJNyF08uz5pGRlUGONYcRESN4buZzDA0fanZZncfWCLs/aRsc/U9wtBqtlnP+Pxh2GQT1MbvCk9Jas99af3yflu/2VlDbNn1oWN9gfjU5nvQEC6P7h8oSxS4m4S6cVmVTJU9teYq1BWsJ8wnjj2l/ZG78XNfY4MvhMO4UzXnLGBzdXAOBfY07RpOvgl5DzK7wpCrqmtmwt8IYJ1dg5VDb9KHoUF/mjuhDWoKxRDHU38VXKzk5CXfhdFodrazZvYZntz5LY2sj1w25jl+P+DUBXs7dkjgtZbtODI6uKQGvABhyodFHj00DN+e7um2y2dlYVElmgXF1vrNt+lCQjwcTB1i4bZqF9EQLMWF+skTRiUi4C6ey8chGMrIzKDhWwPg+43kk9RHiQ7r5drO1RyHvXaPtcng7KHdImAGzfg+DzgMv59oSweHQ7Dxcw/oCY0VLdlHl8elDo2NCWXjOQCYlWEiODnG66UPihA6Hu1LKHdgElGqt53a8JNETHak/wvJNy/m06FP6+vflyalPMiNmRve9EmypN/ZFz3nL2Cdd241b/89dDMMuhYBIsyv8L6VVjcf3afl2bwWV9cYSxYG9Arh2XH/S26YP+XvL9WB30Rn/p+4G8oHudTuccAot9hZe2/kaq3NW49AObh1xK/OHzcfXoxsujXPYYf83bYOj/w4tdRDcD9LuMdouEYPMrvC4miYb/9lbcXxL3P3WegAiA72ZOjCCtERjmWJkUPeaPiRO6FC4K6WigfOBPwL3dUpFosf4puQblmQvobi2mBkxM3hg7ANEBUSZXdaZO5JnXKHnvgu1h8E7GIZd0jY4egK4mf8BsM3uYGtxlXF1Xmhl+8ET04fGxYVx7Xjj6jyxm08fEid09Mp9JfAgcNJRNkqpBcACgJiYmA4+nXAFxTXFLNm4hG9KviE2KJbnZz7PxKiJZpd1ZmoOQe47xmZdR/PAzcMYGJ2cAQPnmD44WmX9+1EAABtOSURBVGtNYdv0oQ2FVr7bV0F92/Sh5OgQbp+WwKQEC6NjQl1q+pA4od3hrpSaC5RprTcrpaae7Pu01quB1QApKSm6vc8nur8GWwN/yf0Lr+54FU83T+4fcz/XJF2DZ3fZsbC5FvI/Mq7S9/0b0BA9Fs5bBkMvAf9wU8srq23i28KK44F+pMZYohgb7sdFo6JIT7QwId5CsF83eb9Fh3Tkyn0SME8pdR7gAwQppf6qtb62c0oTrkJrzadFn7Js0zLKGsq4IP4C7h1zLxF+EWaX9vPsrbDv67bB0f8AWwOExsKUB40+evgA00prbLGTtf9E33zXEWP6UIifJ5MGWI73zfuFOddqHNE12h3uWutHgEcA2q7cF0qwix/bc2wPGVkZbDq6iaSwJJZNWcaoyFFml3VqWhtLFr8fHF1fBj4hMOIqI9D7jTNloy67Q5NXWt12N2g5Ww5U0WJ34OXuRkpsKA+eO4j0hAiG9A2SJYpC1rmLs6O6uZrntj3Hmt1rCPAK4LHxj3Fp4qW4O+FNOsdVHYTct40+evkucPeCgbONQE88Bzy6fihEcUUD6wvL2VBoZUNhBdWNNgCS+gRxw6RYJiVYSI0Nw9fLid9XYYpOCXet9dfA151xLtG9ObSDtQVreWrLU1S3VHP5wMu5Y+QdhPg46QDmpmrj9v+cNVCUCWhjhcvcJ43B0X5hXVpOVUML/9lbwfq2Qc/FlQ0A9A7yYdaQXqQnWpg4wEJEoEwfEqcmV+6i0+SU57AoaxE7KnYwKnIUj457lMFhg80u63/ZbVD4r7bB0Z9AaxOEDYBpjxoDL8LiuqyU5lY7Ww5UGaPkCozpQw4NAd4ejI8P48ZJsaQlRjAgQqYPiTMj4S46zNpo5aktT/FB4QdE+EaQkZ7B+XHnO1cYaQ2lW9oGR78HDRXgFw6jrzPaLlFjuqSPrrVm99Ha4/u0ZO+vpNFmx91NMbJfCHdOb5s+1C8ET5k+JDpAwl20m81h461db/Hctudosjcxf9h8bkm+BX9Pf7NLO+FY0YnB0RWF4O4Ng88zAj1hZpcMjj5a03R8eWJmoZXy2mYA4iP8uSIlmkkJFsYPCCdIpg+JTiThLtol63AWi7MXU1hVyKS+k3go9SHigruunXFKjcdgxwdGoBf/x3gsNh0m3W3swOgTfFafvr65laz9xnrzzAIrBW3Th8L9T0wfmpRoISqkG26xILoNCXdxRg7XHWbppqV8ceALogKieGraU0zrN838FkxrCxR8brRd9nwG9hawDIIZv4XhV0BIv7P31HYHOaXVxnrzAitbio/R6tB4e7iRGhfGZWOiSUu0kNQ7SKYPiS4j4S5OS7O9mVfyXuGF3BfQaG4feTs3DL0BHw8Tb7PXGg5mnxgc3XgM/CNg7M2QfAX0GXlW+uhaa4oqGo7voviffRXUNhnTh4b2DeLm9HjSEy2MkelDwkQS7uKUtNZ8ffBrntj4BCV1JczqP4uFKQvpG9DXvKIq9hotl5w1Rk/dwxeS5hp99PhpZ2VwdGV9i9Ezb7sbtLSqEYCoEF/OH96HSQkWJiVYCJPpQ8JJSLiLkyqqLmLxxsVsKN1AfHA8fznnL4zvM96cYuorYMf7RqCXbAQUxE+BKQ9B0gXgfdK969qlyWZnU9ExMgutZBaWs+NQDVpDoI8HEweE8+sp8aQlRhAbLtOHhHOScBf/o95Wz/M5z/P6ztfxcffhgZQH+EXSL/B06+LVHLYm2POpEegFnxuDoyOHwqw/GOvRgzrvXw8Ohyb/SM3xK/Ps/ZU0tzrwcDOmD907cyBpiRaSo4LxkCWKohuQcBfHaa35eP/HrNi0grLGMi4ccCH3jLkHi6+l64pwOIwVLjlvwY4PobkaAnrD+FuNtkvv4Z32VIeqGo315oVWvi20UtE2fSgxMoCrx8WQnmhhXFy4TB8S3ZL8rhUA7K7czaKsRWwp28KQ8CGsmLaCEREjuq6A8j0nBkdXF4OnPwyZZ3wwGjelUwZH134/fahtvfm+cmP6UESgN5MHRhhLFBMs9A6W6UOi+5Nw7+Gqm6t5ZuszvLPnHYK9gnl8wuNcnHBx12zwVVdu3C2a8xYc2grKDQZMhxmPweDzwatjN0PZ7A62Haw63mrZdrAKu0Pj6+lOalwYV6fGkJZoYVCvQOmbC5cj4d5D2R123it4j2e2PkNNSw1XDrqS20feTrD32b3Bh5YG2P2x0Ucv/JcxOLp3MsxeBMMug8Be7T611pq95fVkFpSTWWjlu32V1DUbSxSTo4KND0ETIhjdPwRvD1miKFybhHsPtK1sG4uyFpFfmc+YXmN4JPURBoWdxeHNDgcUrTcCfec6aKmFoGiYdJfRR49MaveprXXNbCi0Hr+9/3C1MX0oJsyPeSP7kp5gYcKAcEL8ZImi6Fkk3HsQa6OVJzc/ybq964j0i+SJyU9wbuy5Z68lcXTnicHRNaXgFQhDLzQCvX9auwZHN7bYyS6qPB7o+YdrAAj29WRSQjh3Jhi985hwmT4kejYJ9x7A5rDxRv4brNq+imZ7MzcNu4kFyQvw8zwLAVh7xAjznLfgSK4xODphJpzzfzDoPPA8s/1U7A7NjkPG9KHMAiubio4dnz40pn8oD8weRFqChWFRwTJ9SIgfkHB3cd8e+pbF2YvZX72f9Kh0Hkp9iP5B/Tv3SVrqfzA4+mvQDug7GuY8AcMuBf8zW0p5sLLheJhv2GulqsGYPjS4dyDXTehPWqKF1Lgw/Lzkt68QJyN/OlxUaV0pyzYu45/F/6RfYD/+NP1PTOk3pfOewGFvGxy9xgh2Wz2ExED6/cZGXREDT/tU1Y02/rPXaLNkFlo5UGFMH+oV5M2MwW3ThxLCiQyUJYpCnC4JdxfT1NrEy3kv82Lei7gpN+4cdSfXD70eb/dOGMumtdFq+X5wdN0RY/vc5MvbBkePP60+ekurgy3Fx44vUcwpqcKhwd/LnfHx4dwwMZa0BAsJkQGyRFGIdpJwdxFaa74s/pKlm5ZSWlfK7NjZLExZSG//3h0/eXUp5L5jhHrZTnDzbBscfQUkzgbPU19Ra63Zc7SurdVSTtb+ShpajOlDI6KDuWNaAmmJEYzsF4KXh9zaL0RnkHB3Afuq97EkewnfHvqWhJAEXjznRVL7pHbspM21bYOj34L96wEN0alw/nIYesnPDo4uq2k63jfPLLRS9v30IYs/l4429jefINOHhDhrJNy7sbqWOp7PeZ6/7vwrvh6+PJz6MFcOuhIPt3b+b7W3wt4vjUDf9TG0NkJoHEx92NioK3zASQ9taGkla19lW9+8nD1HjelDoX6eTEqwkJ5o3NofHSpLFIXoChLu3ZDWmo/2fcSKzSuwNlq5JPES7hp1F+G+4e05mXHrf84aYyuA+nLwDYVR1xh99OixPznwwu7Q5JScuLV/S/ExbHaNl4cbqbFhXDI6mrQEC0P6yPQhIcwg4d7N7KzYSUZWBtvKtzEsfBhPT3ua4RHt2CmxqvjE4GjrHnD3gkFz2gZHzwKP/76jU2vNgYoTSxS/3WulpqkVMKYP3TgpjrREC2Njw2T6kBBOQMK9m6hqquLprU/z7p53CfUJ5Q8T/8CFCRfips7gA8jGKtj5oRHoBzYYj8VMhAtuhyEXgW/If337sfoWvt1bQWahMU6u5JgxfahvsA/nDutNWmIEkwaEEx7QCStxhBCdSsLdydkddt7d8y7PbHuGupY6rkm6hltH3kqQV9DpnaC1BQr/afTRd38K9mYIT4TpvzH66KGxx7+1udXO5qJjrG+7Os87VG1MH/L2YPyAcBZMjictwUKcxV+WKArh5CTcndiWo1vIyM5gV+UuxvYeyyOpj5AYmvjzB2oNJZtO9NEbK8HPAinzjeWLfUeDUjgcml2Hao5fmW8sqqTJZkwfGhUTwj0zjOlDI6Jl+pAQ3Y2EuxMqayhjxeYV/GPfP+jl14ulU5Yyu//sn79artx/oo9euRc8fIx90ZOvNPZJd/fkcHUj6zeXsKHQ2EXRWmdMH0qIDOCqsTGkJVgYPyCcAJk+JES3Jn+CnYjNbuP1/Nd5fvvz2Bw2fjX8V9w8/OZTb/DVUAk71hqBfjALUBCbBun3QdI86pQf3+2tIPMfe1hfUM7etulDlgAvJiVYSEuwkJZooU/wmW3oJYRwbhLuTmJD6QYWZy+mqKaIqdFTeXDsg/QL6vfT39zaDHs+MwJ9z2fgsEFEEsz8Ha1DLmV7bYCx3vzlPLYdrKLVofHxdCM1Lty4Om+bPiRLFIVwXe0Od6VUP+A1oBeggdVa66c6q7Ce4mDtQZZuXMpXB7+if1B/npvxHOnR6f/7jVpD8XdGoO9YC01VENALnbqAkph5fHmsF+sLK8j6Ip/atulDw6OCjQ9BEy2MjgmVJYpC9CAduXJvBe7XWm9RSgUCm5VSX2itd3ZSbS6tsbWRF3Nf5OW8l3F3c+fu0Xdz3ZDr8HL/0cQga6ER6DlroOoAePrRnDCHzSGzWVeTwDdbqzj0dQVQQb8wX+aO6EtagoWJA8IJ9ZfpQ0L0VO0Od631YeBw29e1Sql8IAqQcD8FrTVfHPiCZZuWcbj+MHPi5nD/mPvp5f+D2aH1Vsh731i+WLoZrdyo6jWBzPgbeLlyGFu2GvubB/lYmTjAwm3TjNv7+4d3bKC0EMJ1dErPXSkVC4wCsn7i1xYACwBiYmI64+m6rb1Ve8nIziDrcBYDQweyKG0RKb1TjF+0NcLuTyBnDbrwnyhHK1b/gXwW/CtWVYyipCgET3fF6JhAFp5jIS0xguEyfUgIcRIdDnelVADwHnCP1rrmx7+utV4NrAZISUnRHX2+7qi2pZZV21fxZv6b+Hr68ui4R7l84OV44GbsuJizBseOD3BrqeWYh4UPHOfzVvMEdjfFMKhXILPHGytaUmPD8JclikKI09ChpFBKeWIE+9+01u93Tkmuw6EdrNu7jpWbV1LZVGls8DX6LsJqy2n6/P+w5byNb8MhGvDhY3sq79vT2OcxkolDenFL2zLFyCCZPiSEOHMdWS2jgBeBfK31is4ryTXssO5gUfYicspzSI5IZmXq7wnduRHbM3Ogfhce2o31juH8Q11Gbf9ZpA7qx+8SLSTK9CEhRCfoyJX7JOCXQK5SalvbY49qrT/ueFndV2VTJU9veZr3C94n2CuEa72mcW7uDoZmXYSHcpDriOPvgQuwJV3M6CGDWBQTKtOHhBCdriOrZTIBucRs0+po5aWcN/hL7nO02BuZWBvE7yr20JvtHFERZPa+BvcRvyB5ZCrD/WT6kBDi7JJP5zqgoaWVrP2VrN35DVmVz9HgUc6ohlYeryyjr92bQ9FzKEu9ht7DptP7NAZHCyFEZ5FwPwN2hya3tJoNhVbWF5RTXLqNUMvbFAdZ6UMrfyirZnz4eAIvXITboDkM+JnB0UIIcbZIuP+M4ooG1heWt00fqsDWWMs57v9hSOQ37Iut5wiKW1r9uWno9fgOvxL82zHqTgghOpmE+49UNRjTh9YXGFviFlc24I6deYG7eTkoi2q/bJaH+nHA05PpvjE8MP43RMdMMrtsIYT4Lz0+3Jtb7Ww+cIwNbdOHckqN6UMB3u5cHnWMyyyZDLJ+zuGWSpb4RfBvn2BifXvx54m/Y1J0mtnlCyHET+px4a61ZvfRWjILrKwvsJK9v5JGmx13N8XIfiE8OimQOY71RBWvQx3aRYO7F8/FJfOKwx9Pd2/uG/Frrk26Fk93WfEihHBePSLcj1Q3kVloJbOgnMzCCqx1zQDER/hzRUo0U2J9mdC0Ht/852BTJqDR/cbxafqvWVaxkaONZcyNn8u9Y+4l0i/S3BcjhBCnwSXDva65lax9Rt88s9BKYVkdAOH+bdOHEi2kxQXT1/ot5CyHv38MrU0QNgCmPUpB/1Qydr/GxpKPGRw2mKVTlzMqcpTJr0oIIU6fS4R7q93B9pJqMts+BN1SfIxWh8bbw43UuDCuSIlmUoKFpF6BuB3ZCtufgS/fgwYr+IbBqF/CiKuoiUjkue2reOubewjwCuCx8Y9xaeKluLvJkAshRPfSLcNda01RRQOZBeWsL7Dyn30V1DYZ04eG9g3i5vR40hMtjOnfNn3o2AHIWQ3vrYGKAnD3hkFzYMRVkDATh5s7HxZ+yMoP7udY0zEuH3g5d466kxCfELNfqhBCtEu3CffK+pbjK1oyC62UVjUCEBXiy/nD+5CWaGHiAAth308fajwG218zJhgV/8d4rH8aTLoLhlwIPsEA5JbnkpGdQa41l5ERI/nzzD+TFJ5kxksUQohO47Th3mSzs6no2PEbiHYcMraKD/TxYOKAcH49dQDpCRb6h/ud2EWxtQXyP2obHP0p2FvAMhCmPwbJV0DIiWEhFY0VPLXlKdYWrsXia2FR2iLmxs+VHRmFEC7BacLd4dDsPFxDZqHRN8/eX0lzqwMPN8Xo/qHcP2sgkxItJEcF4+H+g31atIaD2bD9LdjxvnHF7h8BKTfBiCuhz0j4QWC3Olp5a9dbPLftORpbG7lh6A3cknwLAV4BJrxqIYQ4O0wN99KqRjYUWFlfaOXbQisV9S0ADOwVwNXjYkhPtDAuLvynpw9V7IWct42r9GP7wcMXBp9v9NHjp4H7/x6TfTibjOwMCqsKmdh3Ig+lPkR8cPzZfplCCNHlujTc7Vrz+Y4jbWvOreyz1gMQEejN5IERpLUtU+x1sulDDZWQ954R6CUbAQVxk2HKg5B0AXgH/uRhR+qPsGzTMj4r+oyogChWTlvJ9H7TpQUjhHBZSuuuG2vq3SdR97l+Jb6e7oyLDyMtwUJ6YgQDe51i+pCtyeif57wNBZ+DwwaRQyD5Shh+OQRHnfT5mu3NvLrjVV7IfQGHdnDT8JuYP3Q+Ph6yW6MQovtQSm3WWqecyTFdeuUeEejNm78az+j+IXh7nGLtuMMBB78z+ug7P4CmagjoDeNuMdouvYf/7HP9++C/WZy9mJK6EmbGzGTh2IVEBZz8LwIhhHAlXRruvYN8mDDgFFviWguMQM99G6qKwdPfaLeMuBLipsBp3Ex0oOYAS7KXsL50PfHB8ayetZoJfSd04qsQQgjnZ/5qmbrytj76W3BoKyg34wPR6Y8ZH5B6+Z/WaRpsDazOWc1rO1/Dy92LhSkLuTrpajzdZIMvIUTPY0642xph1z+MPnrhP0HboXcynPNHGH4ZBPY+7VNprflk/ycs37ycsoYy5g2Yx71j7sXiazmLL0AIIZxb14Z7cx18cDvs/BBaaiEoCibeafTRI8/8rtDdlbvJyM5g89HNJIUlsXzKckZGjjwLhQshRPfSteFeUQA764zb/0dcaWwH0I7B0dXN1Ty77VnW7F5DkFcQv53wWy5JuEQ2+BJCiDZdG+6hsbAwF7z82nW43WFnbeFant7yNNUt1Vwx8AruGHUHwd7BnVunEEJ0c10b7r6h7Q727eXbWZS1iJ0VOxkdOZpHxz3KoLBBnVygEEK4BvNXy/wMa6OVJzc/ybq964j0jWRx+mLOiztP7i4VQohTcNpwtzlsvJH/Bn/e/mea7E3cOOxGbkm+BT/P9l35CyFET+KU4f7d4e/IyMpgX/U+0qLSeGjsQ8QGx5pdlhBCdBtOFe6H6g6xbNMyvjjwBdEB0Twz/RmmRE+RFowQQpwhpwj3ptYmXt7xMi/lvgTAnaPu5Pqh1+Pt7m1yZUII0T2ZGu5aa746+BVPbHyC0rpSzul/DgtTFtInoI+ZZQkhRLdnWrjvr97PkuwlbDi0gYSQBF445wXG9RlnVjlCCOFSOhTuSqlzgacAd+AFrfXinzum3lbP89uf5/X81/Fx9+GhsQ9x5eArZYMvIYToRO0Od6WUO/AsMAsoATYqpdZprXee7Jjq5mouWHsB5Y3lXJxwMXePvptw31NsASyEEKJdOnLlngoUaq33ASil3gIuBE4a7iV1JaT4pbBy2kqSI5I78NRCCCFOpSPhHgUc/MHPS4D/aZorpRYACwAi+0fyt/P/hps6883ChBBCnL6znrJa69Va6xStdUo/Sz8JdiGE6AIdSdpSoN8Pfh7d9pgQQgiTdSTcNwKJSqk4pZQXcBWwrnPKEkII0RHt7rlrrVuVUncAn2EshXxJa72j0yoTQgjRbh1a5661/hj4uJNqEUII0Unk000hhHBBEu5CCOGCJNyFEMIFSbgLIYQLUlrrrnsypWqB3V32hM7NAljNLsJJyHtxgrwXJ8h7ccIgrXXgmRzQ1Vv+7tZap3TxczolpdQmeS8M8l6cIO/FCfJenKCU2nSmx0hbRgghXJCEuxBCuKCuDvfVXfx8zkzeixPkvThB3osT5L044Yzfiy79QFUIIUTXkLaMEEK4IAl3IYRwQV0S7kqpc5VSu5VShUqph7viOZ2RUqqfUuorpdROpdQOpdTdZtdkNqWUu1Jqq1LqI7NrMZNSKkQp9a5SapdSKl8pNcHsmsyilLq37c9HnlLqTaWUj9k1dRWl1EtKqTKlVN4PHgtTSn2hlCpo+zH0dM511sP9B4O05wBDgF8opYac7ed1Uq3A/VrrIcB44PYe/F58724g3+winMBTwKda68HACHroe6KUigLuAlK01sMwthO/ytyqutQrwLk/euxh4F9a60TgX20//1ldceV+fJC21roF+H6Qdo+jtT6std7S9nUtxh/gKHOrMo9SKho4H3jB7FrMpJQKBiYDLwJorVu01lXmVmUqD8BXKeUB+AGHTK6ny2itvwEqf/TwhcCrbV+/Clx0OufqinD/qUHaPTbQvqeUigVGAVnmVmKqlcCDgMPsQkwWB5QDL7e1qF5QSvmbXZQZtNalwDKgGDgMVGutPze3KtP10lofbvv6CNDrdA6SD1RNoJQKAN4D7tFa15hdjxmUUnOBMq31ZrNrcQIewGhgldZ6FFDPaf7T29W09ZMvxPgLry/gr5S61tyqnIc21q6f1vr1rgh3GaT9A0opT4xg/5vW+n2z6zHRJGCeUqoIo1U3XSn1V3NLMk0JUKK1/v5fce9ihH1PNBPYr7Uu11rbgPeBiSbXZLajSqk+AG0/lp3OQV0R7jJIu41SSmH0VfO11ivMrsdMWutHtNbRWutYjN8TX2qte+QVmtb6CHBQKTWo7aEZwE4TSzJTMTBeKeXX9udlBj30w+UfWAdc3/b19cCHp3PQWd8VUgZp/5dJwC+BXKXUtrbHHm2bRSt6tjuBv7VdAO0D5ptcjym01llKqXeBLRiry7bSg7YhUEq9CUwFLEqpEuBxYDHwtlLqJuAAcMVpnUu2HxBCCNcjH6gKIYQLknAXQggXJOEuhBAuSMJdCCFckIS7EEK4IAl3IYRwQRLuQgjhgv5/xOQ7d9UTodMAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "geom_solution = deepcopy(PFexample.solution)\n", + "PFexample_hyperbolic = deepcopy(PFexample)\n", + "PFexample_hyperbolic.geometric_solution = geom_solution\n", + "PFexample_hyperbolic.Hyperbolic_beta=0.7\n", + "PFexample_hyperbolic.solve()\n", + "PFexample_hyperbolic.unpackcFunc()\n", + "print(\"Plot of Consumption Functions\")\n", + "for T in range(3):\n", + " x = np.linspace(0,10, 1000,endpoint=True)\n", + " y = PFexample_hyperbolic.cFunc[T+1](x)\n", + " plt.plot(x,y, label = T+1)\n", + "plt.xlim([0, 10])\n", + "plt.legend(title='Period')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "According to `quasi_hyperbolic_algebra`:\n", + "\n", + "$$C_3=X_{3}$$\n", + "\n", + "$$C_{2}=\\frac{1}{1+\\beta\\delta}X_{2}+\\frac{1}{1+\\beta\\delta}$$\n", + "\n", + "$$C_{1}= \\frac{1}{1+\\beta\\delta+\\beta\\delta^2}X_{1}+\\frac{2}{1+\\beta\\delta^2+\\beta\\delta^3}$$\n", + "\n", + "Checking if the slopes match the polynomial expressions:" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1.0\n", + "0.6134969325043818\n", + "0.4551661356268127\n" + ] + } + ], + "source": [ + "print(PFexample_hyperbolic.cFunc[3].derivative(1))\n", + "print(PFexample_hyperbolic.cFunc[2].derivative(1))\n", + "print(PFexample_hyperbolic.cFunc[1].derivative(1))" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1.0\n", + "0.6134969325153374\n", + "0.4551661356395084\n" + ] + } + ], + "source": [ + "print(1/np.polyval([1],0.9))\n", + "print(1/np.polyval([0.7,1],0.9))\n", + "print(1/np.polyval([0.7,0.7,1],0.9))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The same slopes are found." + ] + } + ], + "metadata": { + "jupytext": { + "notebook_metadata_filter": "all" + }, + "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.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/HARK/core.py b/HARK/core.py index 204ff22f6..467616f0a 100644 --- a/HARK/core.py +++ b/HARK/core.py @@ -883,7 +883,11 @@ def solveOneCycle(agent, solution_last): else: T = 1 - solve_dict = {parameter: agent.__dict__[parameter] for parameter in agent.time_inv} + solve_dict = { + parameter: agent.__dict__[parameter] + for parameter in agent.time_inv + if parameter != "geometric_solution" + } solve_dict.update({parameter: None for parameter in agent.time_vary}) # Initialize the solution for this cycle, then iterate on periods @@ -914,7 +918,12 @@ def solveOneCycle(agent, solution_last): # Solve one period, add it to the solution, and move to the next period solution_t = solveOnePeriod(**temp_dict) solution_cycle.insert(0, solution_t) - solution_next = solution_t + # solution_next = s/olution_t + + if hasattr(agent, "geometric_solution"): + solution_next = agent.geometric_solution[T - t - 1] + else: + solution_next = solution_t # Return the list of per-period solutions return solution_cycle diff --git a/examples/Journeys/Quickstart_tutorial/Jounery_1_param.py b/examples/Journeys/Quickstart_tutorial/Jounery_1_param.py index cf012e692..4c072bd69 100644 --- a/examples/Journeys/Quickstart_tutorial/Jounery_1_param.py +++ b/examples/Journeys/Quickstart_tutorial/Jounery_1_param.py @@ -22,6 +22,8 @@ PermGroFacAgg = 1.0 # Aggregate permanent income growth factor (only matters for simulation) T_age = None # Age after which simulated agents are automatically killed T_cycle = 1 # Number of periods in the cycle for this agent type +Hyperbolic_beta = 1 # Naive hyperbolic discount factor + # Make a dictionary to specify a perfect foresight consumer type init_perfect_foresight = { 'CRRA': CRRA, @@ -36,7 +38,8 @@ 'pLvlInitStd' : pLvlInitStd, 'PermGroFacAgg' : PermGroFacAgg, 'T_age' : T_age, - 'T_cycle' : T_cycle + 'T_cycle' : T_cycle, + 'Hyperbolic_beta' : Hyperbolic_beta } # -----------------------------------------------------------------------------