ipynb/kerk_phillips-ak_fixed_point_solver.ipynb

{
    "cells": [
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "# Using the Auberbach-Kotlikoff Fixed-Point Solver with Adaptive Dampening\n",
                "\n",
                "### by [Kerk L. Phillips](https://sites.google.com/site/kerkphillips/home), July 2018\n",
                "\n",
                "The code in this Jupyter notebook was written using Python 3.6."
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "This notebook explains how to use the Auberbach-Kotlikoff fixed-point solver where the dampeing parameter changes as one iterates.  This code and examples can be downloaded at the following Github site.  [https://github.com/kerkphil/DSGE-Utilities](https://github.com/kerkphil/DSGE-Utilities).  Alternatively, this notebook can be found in the public GitHub repository [https://github.com/kerkphil/Jupyter_Notebooks](https://github.com/kerkphil/Jupyter_Notebooks)."
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Auerbach and Kiotlikoff (1987) showed how to solve a deterministic overlapping generations (OLG) model using a fixed-point solver on the time path of aggregate capital for an economy.  Their algorithm proceeds as follows:\n",
                "\n",
                "1) Set the counter $i=0$.\n",
                "\n",
                "2) Take an initial guess for the time path of aggregate capital over a long sample, $\\{K^i_t\\}_{t=1}^T$.  $T$ needs to be large enough that the economy is effectively at the steady state by time $T$.\n",
                "\n",
                "3) Calculate the time paths for wages and interest rates $\\{w_t, r_t\\}_{t=1}^T$ given aggregate capital.\n",
                "\n",
                "4) For all households in all periods, calculate optimal savings given these prices.\n",
                "\n",
                "5) Aggregate the savings of all households each period to get a new measure of the capital stock each period, $\\{\\hat K_t\\}_{t=1}^T$.\n",
                "\n",
                "6) Calculate the difference between $\\{K^i_t\\}_{t=1}^T$ and $\\{\\hat K_t\\}_{t=1}^T$.\n",
                "\n",
                "7) If the difference is close to zero, end the solver, the transition path has been found.  Otherwise update the guess for the time path using $\\{K^{i+1}_t\\}_{t=1}^T = \\xi \\{\\hat K_t\\}_{t=1}^T + (1-\\xi) \\{K^i_t\\}_{t=1}^T$, increment $i = i + 1$ and return to step 2)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "It is well-known that for a value of $0 < \\xi \\le 1$ sufficiently small this proceedure will eventually converge to a fixed-point.\n",
                "\n",
                "However, the smaller $\\xi$ is the longer the algorithm takes to converge.  If $\\xi$ is too large, the algorithm may diverge.\n",
                "\n",
                "One way to speed up convergence is to adopt a value of $xi$ that changes with each iteration.  A simple way to do this is to examine the difference $\\Delta^i \\equiv || \\{K^i_t\\}_{t=1}^T - \\{\\hat K_t\\}_{t=1}^T ||$ each iteration.  If $\\Delta$ rises from one iteration to the next, we will reduce the dampener as this is signal that the algorithm may be diverging.  If the value of $\\Delta$ falls, we can increase the dampener slightly to speed up convergence.\n",
                "\n",
                "$$ \\xi^{i+1} = \\left\\{ \\begin{matrix} \\lambda_1 \\xi^i & \\text{if } \\Delta^i > \\Delta^{i-1} \\\\ \\max(\\lambda_2 \\xi^i, 1) & \\text{if } \\Delta^i \\le \\Delta^{i-1} \\end{matrix} \\right.;\\; \\lambda_1 < 1 < \\lambda_2 $$"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "The following function implements this idea.\n",
                "\n",
                "The code for this is available at [https://github.com/kerkphil/DSGE-Utilities](https://github.com/kerkphil/DSGE-Utilities), but we reproduce that code in the cell below."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 1,
            "metadata": {},
            "outputs": [],
            "source": [
                "import numpy as np\n",
                "\n",
                "def AKsolve(Xguess, funcname, fparams, ccrit, damp, maxiter, shrinkon, \\\n",
                "    shrink, expandon, expand, disttype, display):\n",
                "    '''\n",
                "    This function performs the Auerbach-Kotlikoff contraction mapping on a \n",
                "    function.\n",
                "    \n",
                "    The inputs are:\n",
                "        Xguess:  An initial guess for the fixed point. Can be a scalar or\n",
                "            matrix.\n",
                "        funcname:  Ahe name of the python function.  It must take Xvalue as an\n",
                "            argument with the same dimensions as Xguess, with fparams as \n",
                "            parameters and return a new value for X, Xnew.        \n",
                "        fparams:  A list of parameters used by funcname\n",
                "        ccrit:  The value for distance between Xvalue and Xnew that indicates\n",
                "            convergence to the fixed point\n",
                "        damp:  The weight put on Xnew relative to Xvalue when moving to the\n",
                "            next iteration; Xvalue = damp*Xnew + (1-damp)*Xvalue.\n",
                "        maxiter:  The maximum number of iterations allowed\n",
                "        shrinkon:  If true, the value of damp is scaled down when the distance\n",
                "            between values of X in an iteration increases.\n",
                "        shrink:  The factor by which damp shrinks.\n",
                "        expandon:  If true, the value of damp is scaled up when the distance\n",
                "            between values of X in an iteration does not increase.\n",
                "        expand:  The factor by which damp expands.\n",
                "        disttype:  Indicator variable for the method used to compute distance\n",
                "            between Xvalue and Xnew\n",
                "            1: root mean squared differences  (default)\n",
                "            2: mean absolute deviation\n",
                "            3: maximum absolute deviation\n",
                "        display:  If true, display iterations.\n",
                "    \n",
                "    The outputs are the fixed point, the last iteration's distanceand the\n",
                "        number of iterations performed\n",
                "    '''\n",
                "    # initialize Xvalue\n",
                "    Xvalue = Xguess\n",
                "    # set initial distance measures\n",
                "    dist = 1.0\n",
                "    distold = 2.0\n",
                "    # set counter\n",
                "    count = 0\n",
                "    # begin AK iterations\n",
                "    print('Performing AK contraction mapping')\n",
                "    while dist > ccrit:\n",
                "        if count > maxiter:\n",
                "            break\n",
                "        Xnew = funcname(Xvalue, fparams)\n",
                "        diff = Xnew - Xvalue\n",
                "        if disttype == 2:\n",
                "            dist = np.mean(np.absolute(diff))\n",
                "        elif disttype == 3:\n",
                "            dist = np.amax(np.absolute(diff))\n",
                "        else:\n",
                "            dist = (np.mean(diff**2))**.5\n",
                "        # check if dist is falling, if not lower value of damp\n",
                "        if (dist > distold) and (shrinkon):\n",
                "            # shrink damp and redo with same Xvalue, do not update count\n",
                "            damp = damp * shrink\n",
                "            Xvalue = damp*Xnew + (1-damp)*Xvalue\n",
                "            distold = dist\n",
                "        else:\n",
                "            # update Xvalue and count\n",
                "            count = count + 1\n",
                "            if expandon:\n",
                "                # expand damp if it is < 1.0\n",
                "                if damp < 1.0:\n",
                "                    damp = damp * expand\n",
                "                else:\n",
                "                    damp = 1.0\n",
                "            # take convex combination for new guess\n",
                "            Xvalue = damp*Xnew + (1-damp)*Xvalue\n",
                "            # replace old dist value\n",
                "            distold = dist\n",
                "        # show progress\n",
                "        if display:\n",
                "            print ('count: ', count, 'distance: ', dist, 'damp: ', damp)\n",
                "    \n",
                "    return Xvalue, dist, count"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "## Implementing the Code for a Solow Growth Model's Steady State"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Below is a function that implements the known transition function for the Solow growth model.  Here $\\delta$ is the depreciation rate, $\\gamma$ is the savings rate, and $\\alpha$ is the capital share in GDP.  These values are included in the list `fparams`."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 2,
            "metadata": {},
            "outputs": [],
            "source": [
                "def example(X, fparams):\n",
                "    # An example using the Solow growth model\n",
                "    Xnew = (1-delta)*X + gamma*X**alpha\n",
                "    \n",
                "    return Xnew"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "The following program implements the AK fixed-point solver for this simple model.  We solve for the steady state value of capital $\\bar X$."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 3,
            "metadata": {},
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "Performing AK contraction mapping\n",
                        "Xfixed:     0.495858909362753\n",
                        "Xbar        0.4958419986580077\n",
                        "distance:   9.577551448036736e-07\n",
                        "iterations: 249\n"
                    ]
                }
            ],
            "source": [
                "import numpy as np\n",
                "import scipy.optimize as opt\n",
                "import matplotlib.pyplot as plt\n",
                "\n",
                "# set parameter values for AKsolve     \n",
                "ccrit = 1e-6   # convegence criterion for A-K solution\n",
                "conv = .25      # initital weight on new values in convexifier\n",
                "maxiter = 10000 # maximum number of iterations allowed in A-K solver\n",
                "shrinkon = 1    # binary flag to implement shrinking conv value\n",
                "shrink = .75    # factor to reduce conv if A-K not converging\n",
                "expandon = 1    # binary flag to implement expanding conv value\n",
                "expand = 1.01   # factor to expand conv if A-K is converging\n",
                "disttype = 1    # distance measure formula\n",
                "display = 0     # display iterations on or off\n",
                "\n",
                "# set parameter values for example\n",
                "delta = .08\n",
                "gamma = .05\n",
                "alpha = .33\n",
                "fparams = {delta, gamma, alpha}\n",
                "\n",
                "# Use AK method to find fixed point\n",
                "Xfixed, distance, count = AKsolve(1.0, example, fparams, ccrit, conv, \\\n",
                "    maxiter, shrinkon, shrink, expandon, expand, disttype, display)\n",
                "    \n",
                "# find closed form solution for SS value\n",
                "Xbar = (gamma/delta)**(1/(1-alpha))\n",
                "\n",
                "print('Xfixed:    ', Xfixed)\n",
                "print('Xbar       ', Xbar)\n",
                "print('distance:  ', distance)\n",
                "print('iterations:', count)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Try varying the values of `expand` and `shrink` to see how the number of iterations changes."
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "## Implementing the Code for an OLG Model's Steady State"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Consider a model where a household of age $s$ solves the following dyanamic program.\n",
                "\n",
                "$$ V^s(k_{st}, \\Theta_t) = \\max_{k_{s+1,t+1}, n_{st}} \\frac{c_{st}^{1-\\gamma} - 1}{1-\\gamma} + \\chi_n \\left( 1 - n_{st}^\\upsilon  \\right)^{\\frac{1}{\\upsilon}} + \\beta E \\left\\{ V^{s+1}(k_{s+1,t+1}, \\Theta_{t+1}) \\right\\} $$\n",
                "\n",
                "Here we use an elliptical utility function for leisure as discussed in Evans and Phillips (2018).   This guarantees an interior solution where $0 < n_{st} < 1$\n",
                "\n",
                "The household budget constraint is:\n",
                "$$ c_{st} = w_t n_{st} + (1+r_t-\\delta)k_{st} - k_{s+1,t+1} $$\n",
                "\n",
                "The Euler equations are:\n",
                "$$ c_{st}^{-\\gamma} = \\beta E \\left\\{ c_{s+1,t+1}^{-\\gamma} (1+r_{t+1}-\\delta) \\right\\} $$\n",
                "$$ c_{st}^{-\\gamma} w_t = \\chi_n n_{st}^{\\upsilon-1}(1-n_{st}^{\\upsilon})^{\\frac{1-\\upsilon}{\\upsilon}} $$\n",
                "\n",
                "Suppose also that in the final period of life there is a bequest motive so that households are willing to die with positive assets.  That problem will be:\n",
                "$$ V^S(k_{St}, \\Theta_t) = \\max_{k_{S+1,t+1}, n_{St}} \\frac{c_{st}^{1-\\gamma} - 1}{1-\\gamma} + \\chi_n \\left( 1 - n_{St}^\\upsilon  \\right)^{\\frac{1}{\\upsilon}} + \\chi_b \\frac{k_{S+1,t}^{1-\\gamma}-1}{1-\\gamma} $$\n",
                "\n",
                "This gives a final intertemporal Euler equation.\n",
                "$$ c_{St}^{-\\gamma} = \\chi_b k_{S+1,t+1}^{-\\gamma} $$\n",
                "\n",
                "From the firm's problem we get:\n",
                "$$ Y_t = K_t^\\alpha L_t^{1-\\alpha} $$\n",
                "$$ r_t = \\alpha \\frac{Y_t}{K_t} $$\n",
                "$$ w_t = (1-\\alpha) \\frac{Y_t}{L_t} $$\n",
                "\n",
                "Market clearing gives:\n",
                "$$ K_t = \\sum_{s=1}^S k_{st} $$\n",
                "$$ L_t = \\sum_{s=1}^S n_{st} $$"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "To solve for the steady state we need to solve for the steady-state holdings of capital and labor for each cohort.  Given a set of guesses for $\\bar w$ and $\\bar r$, we can solve for the optimal steady state consumptions and labors by age, $\\{\\bar c_s, \\bar n_s\\}_{s=1}^S$.  The budget constraint then gives us $\\{\\bar k_s\\}_{s=1}^S$, and we can aggregate this and $\\{\\bar n_s\\}_{s=1}^S$ to get aggregate capital and labor.  From there we can update our guesses for $\\bar w$ and $\\bar r$.\n",
                "\n",
                "We need to be able to solve for the value of $\\bar c_1$ that sets $\\bar k_{S+1} = 0$ and in order to do this we have to be able to solve jointly for $\\bar c_s$ and \\bar n_s$ for each age.\n",
                "\n",
                "The function below takes values of $c$ and $n$ and returns tge labor/leisure Euler equation deviation."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 12,
            "metadata": {},
            "outputs": [],
            "source": [
                "# This is the labor-leisure Euler equation set up to be zero by definition\n",
                "def LLEuler(n, c, *HHparams):\n",
                "    # unpack HHparams\n",
                "    wbar = HHparams[0]\n",
                "    \n",
                "    return wbar*c**(-gamma) - chin*n**(ups-1)*(1-n**ups)**((1-ups)/ups)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "We use this function to help in solving for the optimal $\\bar c_1$ in the functions below."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 13,
            "metadata": {},
            "outputs": [],
            "source": [
                "def conshist(c1, *HHparams):\n",
                "    # unpack HHparams\n",
                "    wbar = HHparams[0]\n",
                "    rbar = HHparams[1]\n",
                "    # initialize vectors for consumption, labor and savings\n",
                "    chist = np.zeros(S)\n",
                "    nhist = np.zeros(S)\n",
                "    bhist = np.zeros(S+1)\n",
                "    # set intitial consumption value\n",
                "    chist[0] = c1\n",
                "    for t in range(0, S):\n",
                "        # set up lambda function for fsolve with n as the input\n",
                "        f1 = lambda n:  LLEuler(n, chist[t], *HHparams)\n",
                "        # solve LLEuler for value of n given c\n",
                "        nhist[t] = opt.fsolve(f1, .99)\n",
                "        if printcheck:\n",
                "            # check that LLEuler is close to zero and report\n",
                "            check = f1(nhist[t])\n",
                "            print(\"nhist\", t, \": \", nhist[t], \" check-n: \", check)\n",
                "        # solve for b given c and n\n",
                "        bhist[t+1] = wbar*nhist[t] + (1.+rbar)*bhist[t] - chist[t]\n",
                "        # if not the final period solve for next period's c from interpemporal\n",
                "        # Euler equation\n",
                "        if t< S-1:\n",
                "            chist[t+1] = chist[t]*(bet*(1.+rbar))**(1/gamma)\n",
                "            \n",
                "    return chist, nhist, bhist\n",
                "\n",
                "\n",
                "# parameters    \n",
                "def findc1(c1, *HHparams):\n",
                "    # This function is for fsolve it takes initial consumption and returns \n",
                "    # the final IT Euler error.\n",
                "    chist, nhist, bhist = conshist(c1, *HHparams)\n",
                "    \n",
                "    return chist[S-1] - chib**(-1.0/gamma)*bhist[S]"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Lastly we use these functions to update our values of $\\bar w$ and $\\bar r$."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 14,
            "metadata": {},
            "outputs": [],
            "source": [
                "# find updated wbar and rbar from initial guesses\n",
                "def updatebar(bar, UPparams):\n",
                "    # unpack the bar vector\n",
                "    wbar = bar[0]\n",
                "    rbar = bar[1]\n",
                "    # set parameters to pass to findc1 in fsolve\n",
                "    HHparams = (wbar, rbar, gamma, chin, chib, ups, bet, S, printcheck)\n",
                "    # solve findc1 for value of c1 \n",
                "    c1opt = opt.fsolve(findc1, .1, args=HHparams)\n",
                "    # check that final savings is close to zero and report\n",
                "    # check = findc1(c1opt, *HHparams)\n",
                "    # print \"check-c: \", check, \"wbar & rbar: \", wbar, rbar\n",
                "    # get the consumer's full history\n",
                "    chist, nhist, bhist = conshist(c1opt, *HHparams)\n",
                "    # sum savings and labor to get aggregate capital annd labor inputs\n",
                "    Kbar = np.sum(bhist)\n",
                "    if Kbar < .01:\n",
                "        Kbar = .01\n",
                "    Lbar = np.sum(nhist)\n",
                "    # solve for the implied values of wages and interest rates\n",
                "    wbarnew = (1-alf)*(Kbar/Lbar)**alf\n",
                "    rbarnew = alf*(Lbar/Kbar)**(1-alf) - delta\n",
                "    # put into barnew array\n",
                "    barnew = np.array([wbarnew, rbarnew])\n",
                "    \n",
                "    return barnew"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "We can now proceed to implement the fixed-point solver."
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 15,
            "metadata": {
                "scrolled": true
            },
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "Performing AK contraction mapping\n"
                    ]
                },
                {
                    "name": "stderr",
                    "output_type": "stream",
                    "text": [
                        "C:\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in power\n",
                        "  \n",
                        "C:\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\minpack.py:161: RuntimeWarning: The iteration is not making good progress, as measured by the \n",
                        "  improvement from the last ten iterations.\n",
                        "  warnings.warn(msg, RuntimeWarning)\n"
                    ]
                },
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "count:  1 distance:  0.2171964689195677 damp:  0.2525\n",
                        "count:  2 distance:  0.05182142178578554 damp:  0.255025\n",
                        "count:  2 distance:  0.17733081239204207 damp:  0.19126875\n",
                        "count:  3 distance:  0.017236860082718814 damp:  0.1931814375\n",
                        "count:  3 distance:  0.10166602565831412 damp:  0.14488607812499998\n",
                        "count:  4 distance:  0.029324178992136355 damp:  0.14633493890625\n",
                        "count:  4 distance:  0.043770236021735014 damp:  0.10975120417968749\n",
                        "count:  5 distance:  0.03222890973879605 damp:  0.11084871622148437\n",
                        "count:  6 distance:  0.0291893258194466 damp:  0.11195720338369922\n",
                        "count:  7 distance:  0.02584698506286097 damp:  0.1130767754175362\n",
                        "count:  8 distance:  0.022924448676721547 damp:  0.11420754317171157\n",
                        "count:  9 distance:  0.020298586201267583 damp:  0.11534961860342868\n",
                        "count:  10 distance:  0.017952307977991768 damp:  0.11650311478946297\n",
                        "count:  11 distance:  0.015856831290687982 damp:  0.1176681459373576\n",
                        "count:  12 distance:  0.013987925342548342 damp:  0.11884482739673118\n",
                        "count:  13 distance:  0.01232314378681462 damp:  0.12003327567069849\n",
                        "count:  14 distance:  0.010842102891462693 damp:  0.12123360842740548\n",
                        "count:  15 distance:  0.009526238562741892 damp:  0.12244594451167953\n",
                        "count:  16 distance:  0.008358677672592834 damp:  0.12367040395679633\n",
                        "count:  17 distance:  0.007324098470657398 damp:  0.1249071079963643\n",
                        "count:  18 distance:  0.006408605491831989 damp:  0.12615617907632795\n",
                        "count:  19 distance:  0.005599612299933493 damp:  0.12741774086709123\n",
                        "count:  20 distance:  0.004885732865881902 damp:  0.12869191827576215\n",
                        "count:  21 distance:  0.004256680655257835 damp:  0.12997883745851976\n",
                        "count:  22 distance:  0.0037031750206853904 damp:  0.13127862583310496\n",
                        "count:  23 distance:  0.0032168544078454383 damp:  0.13259141209143602\n",
                        "count:  24 distance:  0.002790195959124682 damp:  0.1339173262123504\n",
                        "count:  25 distance:  0.002416441117127145 damp:  0.1352564994744739\n",
                        "count:  26 distance:  0.0020895268576678065 damp:  0.13660906446921864\n",
                        "count:  27 distance:  0.0018040222050744186 damp:  0.1379751551139108\n",
                        "count:  28 distance:  0.0015550696990500033 damp:  0.13935490666504993\n",
                        "count:  29 distance:  0.001338331501609776 damp:  0.14074845573170044\n",
                        "count:  30 distance:  0.00114993984784601 damp:  0.14215594028901746\n",
                        "count:  31 distance:  0.000986451557177134 damp:  0.14357749969190764\n",
                        "count:  32 distance:  0.0008448063373331728 damp:  0.14501327468882672\n",
                        "count:  33 distance:  0.0007222886256830513 damp:  0.14646340743571498\n",
                        "count:  34 distance:  0.0006164927231406875 damp:  0.14792804151007213\n",
                        "count:  35 distance:  0.0005252909914503322 damp:  0.14940732192517286\n",
                        "count:  36 distance:  0.00044680489134309626 damp:  0.1509013951444246\n",
                        "count:  37 distance:  0.0003793786564347294 damp:  0.15241040909586887\n",
                        "count:  38 distance:  0.0003215554012153618 damp:  0.15393451318682755\n",
                        "count:  39 distance:  0.0002720554811706116 damp:  0.1554738583186958\n",
                        "count:  40 distance:  0.00022975692351706473 damp:  0.15702859690188278\n",
                        "count:  41 distance:  0.00019367776779132188 damp:  0.15859888287090163\n",
                        "count:  42 distance:  0.00016296015457524727 damp:  0.16018487169961065\n",
                        "count:  43 distance:  0.00013685602138535795 damp:  0.16178672041660677\n",
                        "count:  44 distance:  0.00011471426340985696 damp:  0.16340458762077284\n",
                        "count:  45 distance:  9.596923469282903e-05 damp:  0.16503863349698056\n",
                        "count:  46 distance:  8.01304682267989e-05 damp:  0.16668901983195036\n",
                        "count:  47 distance:  6.67735040882687e-05 damp:  0.16835591003026987\n",
                        "count:  48 distance:  5.5531723784212544e-05 damp:  0.17003946913057258\n",
                        "count:  49 distance:  4.608909271863567e-05 damp:  0.1717398638218783\n",
                        "count:  50 distance:  3.817372623065395e-05 damp:  0.1734572624600971\n",
                        "count:  51 distance:  3.155219402129549e-05 damp:  0.17519183508469807\n",
                        "count:  52 distance:  2.6024492070110045e-05 damp:  0.17694375343554505\n",
                        "count:  53 distance:  2.141961031389702e-05 damp:  0.17871319096990052\n",
                        "count:  54 distance:  1.75916363291065e-05 damp:  0.18050032287959952\n",
                        "count:  55 distance:  1.4416335451187893e-05 damp:  0.18230532610839553\n",
                        "count:  56 distance:  1.1788157463247465e-05 damp:  0.18412837936947948\n",
                        "count:  57 distance:  9.617620946635752e-06 damp:  0.1859696631631743\n",
                        "count:  58 distance:  7.829033708107654e-06 damp:  0.18782935979480603\n",
                        "count:  59 distance:  6.3585103954147885e-06 damp:  0.18970765339275408\n",
                        "count:  60 distance:  5.152251612973279e-06 damp:  0.19160472992668162\n",
                        "count:  61 distance:  4.165055444034738e-06 damp:  0.19352077722594843\n",
                        "count:  62 distance:  3.3590303664264097e-06 damp:  0.1954559849982079\n",
                        "count:  63 distance:  2.7024876105670052e-06 damp:  0.19741054484818998\n",
                        "count:  64 distance:  2.1689878970715242e-06 damp:  0.19938465029667188\n",
                        "count:  65 distance:  1.736524999977242e-06 damp:  0.2013784967996386\n",
                        "count:  66 distance:  1.3868260826641255e-06 damp:  0.203392281767635\n",
                        "count:  67 distance:  1.1047563939112341e-06 damp:  0.20542620458531136\n",
                        "count:  68 distance:  8.778103940342298e-07 damp:  0.20748046663116448\n",
                        "distance:   8.778103940342298e-07\n",
                        "iterations: 68\n"
                    ]
                }
            ],
            "source": [
                "# set model parameter values\n",
                "S = 100         # max age\n",
                "alf = .33       # capital share\n",
                "bet = .99       # annual subjective discount factor\n",
                "delta = .08     # annual depreciation rate\n",
                "gamma = 3.0  # intertemporal elasticity of substitution\n",
                "chin = 10.0     # utility weight on leisure\n",
                "ups = 2.2926    # curvature in elliptical utility\n",
                "chib = .1       # utility eight on bequest in last period of life\n",
                "printcheck = 0  # if 1 display each of the labor/leisure Euler deviations\n",
                "\n",
                "# set AKsolve parameter values\n",
                "ccrit = 1e-6   # convegence criterion for A-K solution\n",
                "conv = .25      # initital weight on new values in convexifier\n",
                "maxiter = 200   # maximum number of iterations allowed in A-K solver\n",
                "shrinkon = 1    # binary flag to implement shrinking conv value\n",
                "shrink = .75    # factor to reduce conv if A-K not converging\n",
                "expandon = 1    # binary flag to implement expanding conv value\n",
                "expand = 1.01   # factor to expand conv if A-K is converging\n",
                "disttype = 1    # distance measure formula\n",
                "display = 1     # display iterations on or off\n",
                "\n",
                "# convert beta and delta to per period values\n",
                "bet = bet**(100./S)\n",
                "delta = 1. - (1. - delta)**(100./S)\n",
                "\n",
                "# set paramters to pass to updatebar\n",
                "UPparams = (gamma, chin, chib, ups, bet, S, printcheck)\n",
                "\n",
                "# set initial guess for wbar and rbar\n",
                "wbar = 1\n",
                "rbar = .02\n",
                "bar = np.array([wbar, rbar])\n",
                "\n",
                "# set parameters to pass to findc1 in fsolve\n",
                "HHparams = (wbar, rbar, gamma, chin, chib, ups, bet, S, printcheck)\n",
                "\n",
                "# run the solver\n",
                "bar, distance, count = AKsolve(bar, updatebar, UPparams, ccrit, conv, maxiter, shrinkon, \\\n",
                "    shrink, expandon, expand, disttype, display)\n",
                "print('distance:  ', distance)\n",
                "print('iterations:', count)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Display the results"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 16,
            "metadata": {},
            "outputs": [
                {
                    "name": "stdout",
                    "output_type": "stream",
                    "text": [
                        "wbar:  1.2276001938785823 rbar:  0.016511666209647027\n",
                        "nhist 0 :  0.5882775319264427  check-n:  -4.440892098500626e-14\n",
                        "nhist 1 :  0.5862515532169315  check-n:  -4.796163466380676e-14\n",
                        "nhist 2 :  0.584225892358793  check-n:  -5.1514348342607263e-14\n",
                        "nhist 3 :  0.582200625138234  check-n:  -5.595524044110789e-14\n",
                        "nhist 4 :  0.5801758269435954  check-n:  -6.039613253960852e-14\n",
                        "nhist 5 :  0.5781515727549874  check-n:  -6.483702463810914e-14\n",
                        "nhist 6 :  0.5761279371341351  check-n:  -7.016609515630989e-14\n",
                        "nhist 7 :  0.5741049942144345  check-n:  -7.549516567451064e-14\n",
                        "nhist 8 :  0.5720828176912234  check-n:  -8.171241461241152e-14\n",
                        "nhist 9 :  0.5700614808122659  check-n:  -8.615330671091215e-14\n",
                        "nhist 10 :  0.5680410563684551  check-n:  -9.414691248821327e-14\n",
                        "nhist 11 :  0.5660216166847301  check-n:  -1.021405182655144e-13\n",
                        "nhist 12 :  0.5640032336112137  check-n:  -1.092459456231154e-13\n",
                        "nhist 13 :  0.5619859785145693  check-n:  -1.1546319456101628e-13\n",
                        "nhist 14 :  0.5599699222695758  check-n:  -1.2256862191861728e-13\n",
                        "nhist 15 :  0.5579551352509249  check-n:  -1.3233858453531866e-13\n",
                        "nhist 16 :  0.5559416873252387  check-n:  -1.4210854715202004e-13\n",
                        "nhist 17 :  0.5539296478433092  check-n:  -1.5276668818842154e-13\n",
                        "nhist 18 :  0.5519190856325581  check-n:  -1.616484723854228e-13\n",
                        "nhist 19 :  0.5499100689897203  check-n:  -1.723066134218243e-13\n",
                        "nhist 20 :  0.5479026656737478  check-n:  -1.8740564655672642e-13\n",
                        "nhist 21 :  0.5458969428989346  check-n:  -1.9806378759312793e-13\n",
                        "nhist 22 :  0.5438929673282656  check-n:  -2.1405099914773018e-13\n",
                        "nhist 23 :  0.5418908050669831  check-n:  -2.2737367544323206e-13\n",
                        "nhist 24 :  0.5398905216563761  check-n:  -2.389199948993337e-13\n",
                        "nhist 25 :  0.5378921820677891  check-n:  -2.5490720645393594e-13\n",
                        "nhist 26 :  0.5358958506968492  check-n:  -2.708944180085382e-13\n",
                        "nhist 27 :  0.5339015913579127  check-n:  -2.886579864025407e-13\n",
                        "nhist 28 :  0.5319094672787286  check-n:  -3.055333763768431e-13\n",
                        "nhist 29 :  0.5299195410953186  check-n:  -3.277378368693462e-13\n",
                        "nhist 30 :  0.5279318748470729  check-n:  -3.4638958368304884e-13\n",
                        "nhist 31 :  0.5259465299720603  check-n:  -3.6770586575585185e-13\n",
                        "nhist 32 :  0.5239635673025497  check-n:  -3.8813396940895473e-13\n",
                        "nhist 33 :  0.5219830470607459  check-n:  -4.1389114358025836e-13\n",
                        "nhist 34 :  0.5200050288547325  check-n:  -4.3787196091216174e-13\n",
                        "nhist 35 :  0.5180295716745965  check-n:  0.0\n",
                        "nhist 36 :  0.5160567338889027  check-n:  0.0\n",
                        "nhist 37 :  0.5140865732410883  check-n:  0.0\n",
                        "nhist 38 :  0.5121191468463433  check-n:  8.881784197001252e-16\n",
                        "nhist 39 :  0.5101545111885527  check-n:  -8.881784197001252e-16\n",
                        "nhist 40 :  0.5081927221174652  check-n:  -8.881784197001252e-16\n",
                        "nhist 41 :  0.5062338348460617  check-n:  -8.881784197001252e-16\n",
                        "nhist 42 :  0.5042779039481161  check-n:  -1.7763568394002505e-15\n",
                        "nhist 43 :  0.5023249833559524  check-n:  8.881784197001252e-16\n",
                        "nhist 44 :  0.5003751263583921  check-n:  0.0\n",
                        "nhist 45 :  0.4984283855988888  check-n:  0.0\n",
                        "nhist 46 :  0.4964848130738524  check-n:  8.881784197001252e-16\n",
                        "nhist 47 :  0.49454446013115616  check-n:  -1.7763568394002505e-15\n",
                        "nhist 48 :  0.4926073774688249  check-n:  8.881784197001252e-16\n",
                        "nhist 49 :  0.49067361513390473  check-n:  0.0\n",
                        "nhist 50 :  0.4887432225215076  check-n:  -8.881784197001252e-16\n",
                        "nhist 51 :  0.4868162483740317  check-n:  0.0\n",
                        "nhist 52 :  0.4848927407805539  check-n:  -8.881784197001252e-16\n",
                        "nhist 53 :  0.48297274717639  check-n:  -8.881784197001252e-16\n",
                        "nhist 54 :  0.48105631434282337  check-n:  -8.881784197001252e-16\n",
                        "nhist 55 :  0.47914348840699633  check-n:  0.0\n",
                        "nhist 56 :  0.4772343148419639  check-n:  -8.881784197001252e-16\n",
                        "nhist 57 :  0.4753288384669045  check-n:  -8.881784197001252e-16\n",
                        "nhist 58 :  0.47342710344748845  check-n:  0.0\n",
                        "nhist 59 :  0.4715291532963982  check-n:  -8.881784197001252e-16\n",
                        "nhist 60 :  0.4696350308739988  check-n:  0.0\n",
                        "nhist 61 :  0.46774477838915657  check-n:  0.0\n",
                        "nhist 62 :  0.46585843740020133  check-n:  0.0\n",
                        "nhist 63 :  0.4639760488160321  check-n:  0.0\n",
                        "nhist 64 :  0.46209765289735977  check-n:  0.0\n",
                        "nhist 65 :  0.46022328925808725  check-n:  0.0\n",
                        "nhist 66 :  0.4583529968668218  check-n:  -8.881784197001252e-16\n",
                        "nhist 67 :  0.4564868140485188  check-n:  -8.881784197001252e-16\n",
                        "nhist 68 :  0.4546247784862524  check-n:  -4.440892098500626e-16\n",
                        "nhist 69 :  0.4527669272231107  check-n:  -4.440892098500626e-16\n",
                        "nhist 70 :  0.4509132966642127  check-n:  4.440892098500626e-16\n",
                        "nhist 71 :  0.4490639225788447  check-n:  0.0\n",
                        "nhist 72 :  0.4472188401027118  check-n:  8.881784197001252e-16\n",
                        "nhist 73 :  0.4453780837403032  check-n:  -4.440892098500626e-16\n",
                        "nhist 74 :  0.4435416873673675  check-n:  4.440892098500626e-16\n",
                        "nhist 75 :  0.44170968423349577  check-n:  4.440892098500626e-16\n",
                        "nhist 76 :  0.43988210696480917  check-n:  -8.881784197001252e-16\n",
                        "nhist 77 :  0.4380589875667476  check-n:  8.881784197001252e-16\n",
                        "nhist 78 :  0.436240357426959  check-n:  0.0\n",
                        "nhist 79 :  0.4344262473182829  check-n:  4.440892098500626e-16\n",
                        "nhist 80 :  0.43261668740182946  check-n:  -4.440892098500626e-16\n",
                        "nhist 81 :  0.4308117072301472  check-n:  0.0\n",
                        "nhist 82 :  0.42901133575048045  check-n:  0.0\n",
                        "nhist 83 :  0.42721560130811037  check-n:  0.0\n",
                        "nhist 84 :  0.42542453164978017  check-n:  -4.440892098500626e-16\n",
                        "nhist 85 :  0.423638153927199  check-n:  -8.881784197001252e-16\n",
                        "nhist 86 :  0.4218564947006233  check-n:  -4.440892098500626e-16\n",
                        "nhist 87 :  0.4200795799425132  check-n:  -8.881784197001252e-16\n",
                        "nhist 88 :  0.4183074350412603  check-n:  0.0\n",
                        "nhist 89 :  0.41654008480498494  check-n:  4.440892098500626e-16\n",
                        "nhist 90 :  0.41477755346540124  check-n:  -4.440892098500626e-16\n",
                        "nhist 91 :  0.4130198646817447  check-n:  0.0\n",
                        "nhist 92 :  0.4112670415447634  check-n:  4.440892098500626e-16\n",
                        "nhist 93 :  0.4095191065807674  check-n:  0.0\n",
                        "nhist 94 :  0.40777608175573565  check-n:  0.0\n",
                        "nhist 95 :  0.40603798847947753  check-n:  -4.440892098500626e-16\n",
                        "nhist 96 :  0.4043048476098453  check-n:  4.440892098500626e-16\n",
                        "nhist 97 :  0.4025766794569988  check-n:  -4.440892098500626e-16\n",
                        "nhist 98 :  0.4008535037877148  check-n:  0.0\n",
                        "nhist 99 :  0.3991353398297435  check-n:  -4.440892098500626e-16\n"
                    ]
                },
                {
                    "name": "stderr",
                    "output_type": "stream",
                    "text": [
                        "C:\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:6: RuntimeWarning: invalid value encountered in power\n",
                        "  \n",
                        "C:\\Anaconda3\\lib\\site-packages\\scipy\\optimize\\minpack.py:161: RuntimeWarning: The iteration is not making good progress, as measured by the \n",
                        "  improvement from the last ten iterations.\n",
                        "  warnings.warn(msg, RuntimeWarning)\n"
                    ]
                }
            ],
            "source": [
                "# recover wbar and rbar and print\n",
                "wbar = bar[0]\n",
                "rbar = bar[1]\n",
                "print(\"wbar: \", wbar, \"rbar: \", rbar)\n",
                "HHparams = (wbar, rbar, gamma, chin, chib, ups, bet, S, printcheck)\n",
                "\n",
                "# check the equations have solved\n",
                "c1opt = opt.fsolve(findc1, .1, args=HHparams)\n",
                "printcheck = 1\n",
                "chist, nhist, khist = conshist(c1opt, *HHparams)\n",
                "\n",
                "# drop last period for savings (it is now zero anyway)\n",
                "khist = khist[0:S]"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Plot the steady state values over ages"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": 17,
            "metadata": {},
            "outputs": [
                {
                    "data": {
                        "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAERCAYAAACXT3dwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xd8XOWZ6PHfo95m1HuxXFVdMaazFNNSKMkuAdLYm11270LIpt3N7s3nhk2y2Wy5SeAC2RACIQ2SACFAWAiB0JttsC3Ltoy71WzJljSjXua5f5wzo7FjY2NLmtHM8/189LHmzDsz7/HYz3vO8zZRVYwxxsSHhEhXwBhjzMyxoG+MMXHEgr4xxsQRC/rGGBNHLOgbY0wcsaBvjDFxxIK+MVFARP5JRO6NdD1M7LOgb6KWiNwgImtFpF9EOkTkv0Xk3EjX61SJyAUi0hp+TFW/pap/Fak6mfhhQd9EJRH5AvA94FtAMVAF3A1cFcl6GTPbWdA3UUdEsoGvAzer6qOqOqCqY6r6hKp+WURSReR7ItLu/nxPRFLd114gIq0i8kUROeDeIfxl2Ht/QEQ2i4hfRNpE5Evu8RtF5JUj6qEissD9/ccicrd7t9EvIq+KSIn72T0islVEloe9dreI/KP7WT0icr+IpIlIJvDfQJn7Pv0iUiYit4nIz8Jef6WINItIr4i8ICJ1R7z3l0Rko4j0icgvRSRter4NE2ss6JtodBaQBvzmGM//b+BMYBmwFFgFfDXs+RIgGygHPgPcJSK57nM/Av5GVT1AI/D8+6jXte7nFAAjwOvA2+7jh4HvHFH+48BlwHxgEfBVVR0ArgDaVTXL/WkPf5GILAIeBP4eKASeAp4QkZQj6nI5MBdYAtz4Ps7DxDEL+iYa5QPdqjp+jOc/DnxdVQ+oahfwz8Anw54fc58fU9WngH6gJuy5ehHxqmqPqr79Pur1G1Vdp6rDOA3SsKr+RFUngF8Cy48of6eq7lPVQ8C/ANef4Od8DPidqj6rqmPAfwLpwNlhZe5Q1Xb3vZ/AaQCNOS4L+iYaHQQKRCTpGM+XAXvCHu9xj4Vef0SDMQhkub9/FPgAsEdEXhSRs95HvfaH/T50lMdZhxdn33vU8b0cdn6qGnDfqzysTGfY7+HnZ8x7sqBvotHrwDBw9TGebwfmhD2uco8dl6quUdWrgCLgMeBX7lMDQEawnIiUvM86H03lMep4vKVtDzs/ERH3vdqmoE4mzlnQN1FHVfuA/4OTi79aRDJEJFlErhCRf8fJd39VRApFpMAt+7P3ek8AEUkRkY+LSLabNvEBE+7TG4AGEVnmdoreNgWncrOIVIhIHvBPOCkgcO4Q8t0O66P5FfBBEblYRJKBL+L0Ibw2BXUyce5Yt8/GRJSqfkdE9uN0nP4c8APrcHLjbwNeYKNb/NfAN0/wrT8J3CkiiUAL8An387aJyNeBP+Ckav4R+JtTPI1fAL/HSdf8NlhHVd0qIg8CO9161Ie/SFVbROQTwP/DSemsBz6sqqOnWB9jkGjbRKWgoECrq6sjXQ1jTklTUxNz5szB6/VGuiomTqxbt65bVQuPVy7qrvSrq6tZu3ZtpKthzCmprq7m7rvvZvXq1ZGuiokTIrLn+KWiMOgbY0w8mQgo2/b72bCvlwQRrj298vgvOgUW9I2ZBrt37450FUwUUlXaeofYsK+PDa29rN/Xy6a2PgZHnfEESyuyLegbY8xs1TMwyobWXja29rFhXy8bWnvp7nf641MSE6gv83LtykqWVeawtDKH6vyM47zjqbOgb4wxU2BodILm9j7W7+tlQ2sfG1t72XNwEAARWFCYxQU1RSytyGZJRQ51pV5SkmZ+1LwFfWOMeZ/GJwK07Pez0Q3u6/f1sW2/n4mAMxqyLDuNJRU5XL+qiqUVOTSWe/GkJUe41g4L+sYY8x5Uld0HB93g7qRqmtv7GB4LAJCdnszSyhxW1xWxpCKHpZXZFHmid9FTC/rGGBOms2/YDe69oSt537CzlFNacgKNZdncsGoOSyuzWVaZQ1VeBs5KGbPDCQV9EbkcuB1IBO5V1W8f8fx3gQvdhxlAkarmuM99msllb7+pqg9MRcWNMeZU9QyMsrGtj41hefgD/hEAkhKEmhIPH1xSFsrDLyrOIilxdq9ec9yg704Tvwu4BGgF1ojI46q6OVhGVT8fVv6zuEvMumuOfA1YibPI1Dr3tT1TehbGGHMc/SPjNLX20dQ2GeD3HRoKPT+vMJNzFhSwpCKbpZU51Jd6SUtOjGCNp8eJXOmvArar6k4AEXkIZ8u6zccofz1OoAdnA4ln3TW/EZFncTZ+ePBUKm2MMe9leGyCzR0+Nro5+I1tfezo6ie46kxFbjpLKrL5+BlzWFKRTWN5Nt4o6WidbicS9Ms5fF3wVuCMoxUUkTk4O/kEdyM62mvLj/K6m4CbAKqqqk6gSsYY4xgdD7DNHUnT1NbLBnckzbg7kqbQk8rSimw+vKSMJZXZLCnPJj8rNcK1jpwTCfpH66E41ipt1wEPuzsJnfBrVfUe4B6AlStXRtcKcMaYqDE+EWBH18BkJ2tbH1s6fIyOT46kWVKRzd/UznNG0lTkUOxNnVUdrdPtRIJ+K4dvBlHBsTesuA64+YjXXnDEa1848eoZY+JVIKDs7B5gU5uzZEFTax/N7T6GxpxryqzUJBrLvdx4djWLy7NZWpFDZV66BfjjOJGgvwZYKCJzcXbuuQ644chCIlID5OLsehT0DPCtsE2pL8VZp9wYY0JUlT0HB9nY1kdTa3AsvI/+kcOHSn7s9EqWVmazuDyHeQWZJCRYgH+/jhv0VXVcRG7BCeCJwH2q2uxuOLFWVR93i14PPKRhC/Sr6iER+QZOwwHOZtWHpvYUjDGziarS2jPkpmecK/hNbX2hsfApSQnUlXr5yIpyGsuzWVKRzYLC2T9UMlpE3SYqK1euVFtP35jYEFxVclNbn9vR6vz0Do4BkJwo1JZ4WVyRzWI3wC8q9pBsAf59E5F1qrryeOVsRq4xZkqoKu19wzS19rrB3UdTay89boAPTna6vKGExRXZLCnPYVFJFqlJsTcWPppZ0DfGvG/hV/DBAL+prY9DA86ywUkJwsJiD5fWl9DoXsXXlnhicrLTbGNB3xjzno6Womlu94UCfGKCsKjYw+q6IhaXZ7O4IscCfBSzoG+MCQl2sk5ewTudrMEUzZEBvrE8m7oYXa4gVlnQNyZOqSp7Dw2yqc0XCu6b2ic7WZPcAG8pmthiQd+YOBAIKLsOOhOdnB8fm9r78LvDJIOjaK5oLKGx3AnwNSUe62SNQRb0jYkxEwFlR1d/KEXT3Oajub2PAXfz7ZSkBOpKPFy5tCwU4BcVeyKydZ+ZeRb0jZnFRscDvHvAT7N75b6prY8tHf7QUgVpyQnUl3r56GkVNJY5OfiFxVk2Dj6OWdA3ZpYYHptga6efTW3Odn3N7T62dvgZnXAWG8tMSaShLJvrVlWGOlnnFWTaTFZzGAv6xkSh/pFxtnT4Qvn35vY+3j3QH9p4Ozs9mYYyL395TjUN5dk0lHmZm29r0Zjjs6BvTIT1DIzS3D6Zntnc7mPXwYHQhh8FWaksLvdySX0xDWVeGsqyqci11STNybGgb8wMUVX2+0Zobp+8em9u99HWO7llX3lOOg1lXq5aVs7iCi+NZdkUedMiWGsTa6ZkY3S3zLXAbTibpGxQ1Rvc4xNAk1tsr6peOQX1NiaqBQLK7oMDNLf73B/nCv6gO4tVBObmZ7JiTi6fOmsOjeXZ1Jd6yc1MiXDNTaybko3RRWQhzjr556hqj4gUhb3FkKoum+J6GxM1QiNo2n1sDgvwwSGSyYnCwiIPF9UW0VDmDc1izUy1G20z86ZqY/S/Bu5S1R4AVT0w1RU1JhoMjIyztdN3WHrm3f39oRE0GSmJ1LlDJIP594XFtpKkiR5TtTH6IgAReRUnBXSbqj7tPpcmImuBceDbqvrYkR9gG6ObaNTdP/InV+/hHax5mSnOCJpzq2koc9IzcwsySbQRNCaKTdXG6EnAQpz9cCuAl0WkUVV7gSpVbReRecDzItKkqjsOezPbGN1EUCCg7OsZdIP75BX8Af9IqEx4B2tDmZeGci8l3jQbQWNmnanaGL0VeENVx4BdItKC0wisUdV2AFXdKSIvAMuBHRgTAaPjAbbt97O5w7mC39zuY0uHD7+7F2tigrCgMItzFxRQX+alvsxLQ2k22RnJEa65MVNjqjZGfwxnj9wfi0gBTrpnp7sh+qCqjrjHzwH+fcpqb8x76Bscc4J7MMB3+Nh+wM/YhHMzGcy/X73cuXqvL/OyqNhWkTSxbao2Rn8GuFRENgMTwJdV9aCInA38QEQCQAJOTn/zMT7KmJMS3OQjGNiDaZrw8e+FnlTqS71cWFNIfZmXulIv1fmWfzfxxzZGN7NKcHjklg6/G+SdDlafu0SwCMwtyKS+1LlyD/5Z5LEJTia22cboZtbrGRhlSzA90+FjS4f/sPRMenIitaUePrS0LBTca0s8ZKTYP2tjjsX+d5iICwScHZy2dPgmg3y7j/a+4VCZIk8qdaVeLqgppL7USc/Y8Ehj3j8L+mZGDY1OsLUzeOXuXL1v7ZicvZqYIMwryOT0uXnUucG9vtRLoSc1wjU3JjZY0DfTQlXp9A2HAnswyO/qnpzc5ElNoq7Uy1+srKSu1ENdqY2eMWa6WdA3p2xkfIJ39/eHAvyWDh9bOn2hDbYBKvPSqS/1cuXSstDVuy0PbMzMs6Bv3pcD/uFQYN/qBvkdXf2Mu5t7pCUnUONusB1Mz9SUePCm2eQmY6KBBX1zVKPjAXZ0OVfvWzv9oU7W7v7RUJnS7DTqSr2sri+irtRLbYl1rhoT7SzoG7r8I25wn0zP7OjqDw2NTElKoKbYw4U1RdSWep38e4mt/W7MbGRBP46MjE+w/UA/Wzv8bO2cvIIPv3ov8aZRW+rhgpoi6ko9oZUjbXNtY2KDBf0YFNyWb0unbzLAH5F7T01KYJF79V5X6qXWrt6NiQsW9Ge5odEJtu33h1IzwSv48JEz5Tnp1JR4WF1fRG2Jk56pzrerd2PikQX9WSIQUFp7hg6/eu/0sztsU4+MlEQWFXu4orGU2hIPNSXO1bstC2yMCbKgH4V6B0fZ2umnpXMyuLd0+hl0Z62KwJy8DGpLvFy1rIzaEmdiU2VuBgk2csYY8x5OKOiLyOXA7ThLK9+rqt8+SplrgdtwdtXaoKo3uMc/DXzVLfZNVX1gCuodE4LDIp3g7gT4lk4/HWFrzuRkJFNb4uHalZXOlXupl0XFWbaomDHmpBw3cohIInAXcAnODllrROTx8HXxRWQh8I/AOaraIyJF7vE84GvASpzGYJ372p6pP5XoFVzvPRjcW9yf8I7V5ERhQZGHM+flT6ZmSr0UeVJt1qoxZsqcyOXiKmC7qu4EEJGHgKuA8M1Q/hq4KxjMVfWAe/wy4FlVPeS+9lngcuDBqal+9OkbGnOD+mRapmW/H7+73js4Hau1JR4uriuipsRDbYmXeYWZJFvHqjFmmp1I0C8H9oU9bgXOOKLMIgAReRUnBXSbqj59jNeWH/kBInITcBNAVVXVidY9ooJj3p2RM5NX7+GpGU9aErUlHq5eVu4Gdw+LbEkCY0wEnUjQP1pu4cjttpJwNkK/AGfj9JdFpPEEX4uq3gPcA87OWSdQpxkTCCj7egYPu2pv6fSzq3uAibDUzPzCLM6Ym0dNiTeUninNTrPUjDEmqpxI0G8FKsMeVwDtRynzhqqOAbtEpAWnEWjFaQjCX/vCyVZ2OqkqXf0joSv2lk4/2/b72ba/n6GxiVC5qrwMako8XN5QwiL36n1ugaVmjDGzw4kE/TXAQhGZC7QB1wE3HFHmMeB64MciUoCT7tkJ7AC+JSK5brlLcTp8I8o/PMa2/X5aOvtDE5u27e/n0MDkcgT5mSnUlHi4blUlNcXOlfuiYg+ZqTZqxhgzex03gqnquIjcAjyDk6+/T1WbReTrwFpVfdx97lIR2QxMAF9W1YMAIvINnIYD4OvBTt2ZMDw2wY6uybz7tk7nyr2tdyhUJjMlkUUlHi6tL2ZR8WTevSDLdmoyxsQeUY2qFDorV67UtWvXvu/XDY6O80JLVygt07Lfz+7uAdy0eyjvHrxir3X/LM9JtwlNxphZT0TWqerK45WLmVzF0OgEf/fztxGB6vxMFhZl8aHFpSwq8VBT7KHa8u7GGBM7QT8/K5UnP3su8wuzSE+xPVaNMeZoYiboAzSWZ0e6CsYYE9WiLqcvIl3AnkjXwxhjZpk5qlp4vEJRF/SNMcZMH+vZNMaYOGJB35jjEJEXROSvIl0PY6aCBX0TV0Rkt4isjnQ9jIkUC/rGzBARianRcmZ2sqBv4p6I5IrIkyLSJSI97u8VRxSbLyJviUifiPzW3SAo+PorRaRZRHrdVFBd2HO7ReQfRGQjMGCB30SaBX1jnP8H9wNzgFeBc4DN7j4PACXAl3CWCn8UaADuEJFCEXkG+A3OnJcrgaeAJ0QkJez9rwc+COSo6jjGRJAFfRP3VPWgqj6iqoPAp4ALcRYOvFVEyoFqnC1DzwYWAG8C1wJ3AHuBR3C2E/0B8J9Auls26A5V3aeqQxgTYXaraeKeiGQA38XZyrME5/9FgvvzSaAXeFdVx0Tk10A9kAysxrn6DwCPA14gE2e3uPAd4sJ3jzMmouxK3xj4IlADfB5nGfBz3OPrgRb39/CNhDzAGE7Avwt4QVWXqWo50O+WbQsrbzMgTdSwK30Tj5JFJC3scS4whHMR1A/8g3v8DCADyAE+KSI/B/4CKAUedsukAx8UkYuBPpyd4kaA16b5HIw5KRb0TTx66ojH9+ME7x/jXL33uMffxLli3wMU4nTyprmPb3XL3AV0hb3nG8CHVXVyGzZjokjUrb1TUFCg1dXVka6GMcbMKuvWres+kQXXou5Kv7q6mpPZOcsYY+KZiJzQ6sRRF/RP1uh4gI/d8zr1pV6WVGSzuDyHhcVZtluWMcaEiZmg3zM4SlpSIo+vb+fnb+4FIDUpgfoyL4vLs1lcns2SihwWFGWRaHviGmPiVNTl9E92Y/SgQEDZc2iQja29NLX2sbGtj+a2PgZGJwBIT06koczL4ors0B3BvIJM2xzdGDOrnejG6DEX9I9mIqDs6u6nqa2PDfv62NTWR3O7j6ExpyHITEmkoTybJeXZbmOQw5y8DGsIjDGzxokG/ZhJ77yXxARhQZGHBUUerlnurKM1EVC2H3AagqbWXja09vHTN/YwMh4AwJOWRGOZezdQkc2S8hwq89IRsYbAGDN7xcWV/okamwjw7v5+mtp63cagjy0dfkYnnIYgOz3Z6R+oyA71E1TkWkNgjIk8S+9MkdHxANv2+2lq62Njax9Nbb20dPoZm3D+3nIzkllckcPici+Ly3NYUpFNaXaaNQTGmBllQX8ajYxP0NLpZ0NrH5ta+2hq66Nlv5+JgPN3WZCVQqPbR9Dojhoq9qZaQ2CMmTaW059GqUmJLKnIYUlFTujY8NgEWzp8bGrrcxqDtj5efrc7rCFIZUlFdqgxWFKRTZE37VgfYYwx08KC/hRJS05keVUuy6ty+aR7bGh0gs0dfaGho02tfbzQcgC3HaDYm+r2DeSwuMJJDxV6UiN2DsaY2GdBfxqlpyRy2pw8TpsT2lmPwdFxNrf72OjeDWxs6+O5rQcIZtlKvGluSig71GlckGUNgTFmaljQn2EZKUmsrM5jZfVkQ9A/EmwIetnU5vQRPLd1f6ghKM1OC6WFgiOH8q0hMMacBAv6USArNYlVc/NYNXeyIfAPj7GpzRdqBJra+nh28/7Q8+U56TSWe1lSkUOjO3w0LzPlaG9vjDEhFvSjlCctmbPm53PW/PzQMd/wGM1tPncegdMgPNNsDYEx5sRZ0J9FvEdpCPqGxmhud/sH3OGjRzYER04oy7WGwJi4ZUF/lstOT+bs+QWcPb8gdKxvaIxmt5O4qc1pEJ5u7gw9H94Q2B2BMfHFgn4Myk5P5uwFBZy94E8bgia3MThaQ9BY7ixDHWwIrLPYmNhjQT9OHK8hCN4RHK2PwBoCY2KHBf04dsyGwO0jaGrz0dTae1hDUOYOHw02BI3l2TahzJhZxIK+Ocyx+gg2t08OH93U1sfvw4aPBieUBe8KFpfbEhPGRCsL+ua4stP/dNSQf3iMZrchONqEsiJPauhOINgQ2KJzxkSeBX1zUjxpyZw5L58z5002BMGZxU3uFpVNbYevNVSQlRq6G2goc+4MynNsPwJjZpIFfTNljjazeHB0nC0dPppa+9jk3hmErz6am5EcuiNoLHPuCGyHMmOmjwV9M60yUpL+ZNG58GWoN7X52NTexw9f2sm42xB405JoKHPmETSUOXcG1fm2eb0xU8GCvplx4ctQB42MT7Ctsz80fHRzex8/fm03o+6exZkpiW5KyEkLNZZnM68gk6TEhEidhjGz0ikFfRG5HLgdSATuVdVvH/H8jcB/AG3uoTtV9d5T+UwTm1KTEp2lIiqyQ8fGJpytKpvdu4Gmtj5+8dYehsechiAtOYG6Ui+Nbv9AQ1k2i4o9pCRZQ2DMsZz0dokikghsAy4BWoE1wPWqujmszI3ASlW95UTfdzZsl2giZ3wiwK7uAXfoqNMYbG730T8yDkBKYgI1JZ5QI9BYnk1tiYe05MQI19yY6TUT2yWuArar6k73Ax8CrgI2v+erjDkFSYkJLCz2sLDYw0dWOMcCAWXPocHQ8NHmdh9PNXXy4Fv7AEhMEBYUZtFQHrwryKa+zEtWqmU3Tfw5lX/15cC+sMetwBlHKfdRETkf567g86q678gCInITcBNAVVXVKVTJxKOEBGFuQSZzCzL58NIyAFSV1p4hmtt9NLupoZff7ebRt9tCr5tbkElDmdM/0FDmNAi2AqmJdacS9I82lOLIXNETwIOqOiIifws8AFz0Jy9SvQe4B5z0zinUyRgARITKvAwq8zK4vLEkdPyAb5hN7X2hfoJ39vby5MaO0PPlOenUuw1AsEGwSWUmlpxK0G8FKsMeVwDt4QVU9WDYwx8C/3YKn2fMKSvypnGRN42LaotDx3oHRydnF7t3Bn/YMjm7uCArhQa3EQhOKqvKy7CGwMxKpxL01wALRWQuzuic64AbwguISKmqBi+jrgS2nMLnGTMtcjJSOGdBAeeELTw3MOJMKgtvDF4Nm0vgSU2irswbaggayrwsKMoi2YaQmih30kFfVcdF5BbgGZwhm/eparOIfB1Yq6qPA7eKyJXAOHAIuHEK6mzMtMtM/dMN7INzCTa199Hc7nQYP/jW3tAQ0pSkBGpLPIc1BLUlXtJTbOSQiR4nPWRzutiQTTObTASUXd39bGrzhRqC5nYffUNjACQIzCvMchuCycYgJ8M6jM3UOtEhmxb0jZliwZFDm9300Ga3MejoGw6VCXYYBxuC+jIvZdlp1k9gTtpMjNM3xhxF+MihyxomRw4d7B8JNQTNR+kwzs1IdhsC526gvtTLvMIsEm3NITOFLOgbM0Pys1I5b2Eh5y0sDB1zViH1h+4GNnf4DltzKC05gZoSb6gRqC/zUmf9BOYUWHrHmCgzNhFgR1c/m9uD6SHnrsA37Cw1kSDBiWVOWijYGBTY/sVxzdI7xsxSyYkJ1JY4I3+CS02oKm29Q5MNQYePdXt6eHzD5NSYYm9qqAGoL3UahDl5GbYktTmMBX1jZgERoSI3g4rcDC4N6yfoHRwNNQRbOpzG4KWwTWoyUxKpLT08PbSo2Bagi2eW3jEmxgyPTbD9QH8oLbSlw8/mjsmVSBME5hdmHZYaqiu19NBsZ+kdY+JUWnJiaAvK4EopgYCyr2fQuRtw00Nrd/fw2/WT6aEiTyp1pU4D4DQIHuYW2OihWGNB35g4kJAgzMnPZE5+Jpc3loaOB9NDmzt8bOnw09zex6vbu0PLTaQlJ1BT7AndDdSXeqkttWWpZzNL7xhjDjMy7qSHtnT4D7szCM4yBpiTn0Fdide9M/BQV+qlItc2tI8kS+8YY05KalKiO0FscutKVaWjb5gtHT73x+kneGZzZ2hymSctyW0IPKE0UY3tWhZ1LOgbY45LRCjLSacsJ52L6yaXpR4YGadlv/+wxuDX61oZHJ0AJucUBBuBYINQ4rUlJyLFgr4x5qRlpiaxoiqXFVW5oWOBgLL3kNNpvKXTaRDW7zt8s5qcjGRqS9w7AjdNtLA4y+4KZoAFfWPMlEpIEKoLMqkuyOSKxZOdxr7hMba6/QRbO527gofe2sfQmHNXkJggzCvIpLbU6zYIdlcwHSzoG2NmhDctmVVz81g1d3KPgonwuwI3PfT2nh6eCJtpHLwrqHX7C2pLnAlmtv7QybGgb4yJmMSwTe0/cMRdQUvn4X0Fv1wzeVcgAnPzM6ktDTYGzt2BjSA6Pgv6xpio401L5vTqPE4P27ksEFD2HBqkxU0Nbe10lp94qqkzVCYrNYmaEo97Z+Ch1h1B5E1LjsRpRCUbp2+MmdUGRsbZ2umnpdNpCLa6DUJwVVJwNq2pLfE4DYJ7VzC3IDOm9jS2cfrGmLiQmZrEaXNyOW3O5Aii4LyCrZ2+UIOwpcPHi9u6QrONUxITmF+URZ3bGNS4o4mKPKkxnSKyoG+MiTnh8wouqp2cVzA67uxVMHlH4Oe1HQd59J22UJmcjGRqip300CK3A7mmxBMzS0/ExlkYY8wJSElKCE0UY/nk8d7B0bAUkZ+WTh+PvN0WWpkUoCJ3MkVUUzJ7U0QW9I0xcS8nI4Uz5+Vz5rz80LHgBvctnf7QrONt+/280DKZIkpOFOYXZrGo2G0M3D+jeRTRKQV9EbkcuB1IBO5V1W8fo9yfA78GTldV66U1xkS98A3uV9dPpohGxifY2TUQagxaOv1/sotZVmoSC4uznBRRsdMYLCrxRMWeBScd9EUkEbgLuARoBdaIyOOquvmIch7gVuDNU6moMcZEg9SkxMkUURj/8Bjb9vtp6eynxe1AfnpTJw++tS9UpiArhUXFntCdQfDPmewvOJVPWgVsV9WdACLyEHAVsPmIct8A/h340il8ljHGRDVPWjKnzcnjtDmTcwtUla7+Eef6scolAAAaZElEQVSuoNPvNgp+frV2X2hROnCGlNaUeFhZncvfXbBgWut5KkG/HNgX9rgVOCO8gIgsBypV9UkROWbQF5GbgJsAqqqqTqFKxhgTPUSEIk8aRZ40zltYGDoeCDgb3W91G4JgY7Bml8IF01unUwn6R+ulCM30EpEE4LvAjcd7I1W9B7gHnMlZp1AnY4yJegkJk/0Fl4T1F8zEZNlTGWvUSnADTkcF0B722AM0Ai+IyG7gTOBxETnujDFjjIlHMzHi56SXYRCRJGAbcDHQBqwBblDV5mOUfwH40vFG74hIF7DnpCpljDHxa46qFh6v0Emnd1R1XERuAZ7BGbJ5n6o2i8jXgbWq+vhJvu9xK22MMebkRN2Ca8YYY6bP7Jo/bEyUEJGPi8jvT+H1L4jIX01lnYw5ERb0TcwTkRtEZK2I9ItIh4j8t4iceyrvqao/V9VLwz5DRWR6B1gbMwUs6JuYJiJfAL4HfAsoBqqAu3EmEhoTdyzom5glItnA14GbVfVRVR1Q1TFVfUJVvywiq0TkdRHpde8A7hSRlLDXq4jcKiI7RaRbRP7DnX+CiNwoIq+4v7/kvmSDezfxMRHJFZEnRaRLRHrc3ytm/C/BmCNY0Dex7CwgDfjNMZ6fAD4PFLhlLwbWicg6EQkOPb4G507BB9yMM+/kTve5JBF5BEh3H/+lqmap6i9x/m/dD8zBubsYAoKvMyZiLOibWJYPdKvq+NGeVNV1qvqGqo6r6m7gB8AeVT0NCE4ivBf4gvv4fwFLwt5iHvBdVT3dffytsPc+qKqPqOqgqvqBfwH+bArPzZiTYuvpm1h2ECgQkaSjBX4RWQR8ByegZ+D8f+gSkQ1hxU4DXlTVQyKyCxgOey4HuDNsFmWWiHhU1S8iGTjLkFwOBPfx84hIoqpOYEyE2JW+iWWv4wTpq4/x/PeBrcBCVfUC9+EE8rNUdalbpj+sfBXQF/ZY3LLL3Mfnulf1AF8EaoAz3Pc+P+w1xkSMBX0Ts1S1D/g/wF0icrWIZIhIsohcISL/jrM+lA/oF5FanPz9uKoOuo8BPgJcKCKNwOc4PGj3ALe4v+8HVoc958HJ4/eKSB7wtWk6TWPeFwv6Jqap6ndwcvJfBbpwlgO/BXgMZ4+HGwA/8EPgx4CIyEacfSAAnsfpqN0AZLqPg1f7O4GVbnkF/q87EuhanM7fdKAbeAN4elpP1JgTFHXLMBQUFGh1dXWkq2EM69ato6GhgbS0tEhXxZjjWrduXfe0Lrg2Xaqrq1m71rbRNZEnIjz22GMsWGATbU30E5ETWp046oK+MbOFquIbGqd7YIRu/wgHB0bpGRyld3CMvqExfENj+IfH8Y+MMzAyzuDoBEOj4wyPBRidCDA6HmBsIsBEQJlQJfymO0EgKSGBxAQhOVFITU4kLTmB9OREMlKSyExNJDMlCW96Mt60ZHIyksnNTCE/M4W8zBQKPakUeVLJSk2akTXazexhQd+YY/APj7Hv0CDPbt5Pa88gHX3DtPcO0dk3zH7/MPt9I4yOB4762tSkBLzpyXjSkvCkJpGZmkRuRgoZKU7wTklKICUxkeQkIVGExASZDM6qBBTGA0pAldHxACPjE4yMBRgam2BgdILBkXG6/YP4h50GZmD06KNAM1ISKclOoyw7nbKcNCpyM6hyd2yaW5BJbkayNQpxxoK+iWvDYxPs7BpgZ3c/O7sG2NU9wO6DA+w5OMihgdHDyqYkJVCWnUZJdhqnVeVS7E2j0JNKoSeV/MxU8tyr7JyMZNKSE2f0PEbHA/QOjnJwYJSD/aN09Q9zwDfCft8Inb4h2nuH+WNLF13+kcNel52ezNyCTBYUZbGoOIuFRR5qSz2UeNOsMYhRFvRNXBifCLD74ABbOvxs7fTR0tnP9gN+9h4aJBCWVinLTqO6IJPLGoqpysukMi+dytwMynPTyc9MidpAmJKUQJE3jSLve3c6D49N0NozyJ6Dg+w+OMgut7F7cVsXD69rDZXLyUimrsRLY7mXxRU5LCnPZk5+RtSevzlxFvRNzBkZn6Cl08/G1j6a2/vY3O5ja6efETcVk5QgzC3IpKEsm6uWlbOgKIv5hVnMLcgkPWVmr9BnWlpyIguKPCwo8vzJcz0Do7x7oJ+tnT42t/vY0uHjgdf3MDq+C3AagmWVOayoymXlnFyWV+XG/N9XLJqxoC8iicBaoE1VPzRTn2tim6qyq3uAd/b2sn6f87O108fYhHP5np2eTEOZl0+dNYe6Ui+1JV7mF2WSmmTB6ki5mSmsmpvHqrl5oWNjEwG27Xca0PV7e3lnXw8vbutC1Wk8F1dkc8bcfM6en8/p1XnWCMwCMzZO313XfCXgfa+gv3LlSrUhm+ZYhscm2Njax5rdh3h7Tw9v7+2hZ3AMgMyURJZU5LC0MoclFdksLs+mIjfdUhJTrG9ojLf39vDWrkO8tesQG/b1Mh5QkhOFFVW5/FlNIecvLKS+1EtCgv3dzxQRWaeqK49bbiaCvruO+AM4Kw1+wYK+OVFDoxO8vbeHN3Ye5M2dh1jf2hsaMTO/MJPT5uSyospJNSwoyiLRgsyMGxgZZ83uQ7y24yCvvNvN5g4fAIWeVC6qKeKiuiLOW1hARoplk6dTtAX9h4F/xVmP5EtHBn0RuQm4CaCqquq0PXtOaI6BiUFjEwHW7+vlte0HeXVHN+/s7WFsQklMEBrLvKyam8fp1XmsrM4jLzPl+G9oZtwB/zAvb+vmjy0HeHFbF/7hcVKSEjh/YSGXNRSzuq6YXPvuplzUBH0R+RDwAVX9OxG5gKME/XB2pR9fVJXdBwd5aVsXL7/bxes7DjIwOoEINJZlc/b8fM6cn8/KObl40pIjXV3zPo1NBFiz6xC/37yf3zd30t43TGKCcM6CAj60pJTL6kvIzrDvdSpEU9D/V+CTwDjOLkZe4FFV/cTRylvQj33DYxO8tqObF1q6eKGli72HBgGoysvgvIUFnLewgDPn5ZOTYVeDsURVaWrr46mmTp7c2E5rzxDJicKFNUV8ZEU5F9YWWQf7KYiaoH/Yh9mVftzq6BviuS0HeG7Lfl7bcZCR8QDpyYmcPT+fC2oKOX9RIXPyMyNdTTNDVJUNrX08vr6dxze0090/gjctiauWlfMXKytYXJ5tHfDv04kGfetZMdNCVdna6ef3zft5dksnm9qczr2qvAyuX1XFRbVFnDEvz67s4pSIsKwyh2WVOfzTB2p5bcdBHnm7lV+t3cdP39hDbYmH606v5JoVFWSnW/pnKkXd0sp2pT97BQLKO/t6eHpTJ88072fvoUFEYEVVLqvrirmkvoj5hVl2BWeOqW9ojCc3tvPLNfvY2NpHWnICH1pSxqfOmsOSipxIVy+qRWV650RY0J9dJgLK2t2HeKqpg6ebO9nvGyE50emou6yhhIvriijy2Hr05v3b1NbHz9/cy2/XtzE4OsHyqhxuPLuaKxpLSUmy/Z+OZEHfTJtAQHl7bw9Pbuzgd00ddPlHSE1K4IKaQq5oLOWiuiK8NtLGTBHf8BiPrGvlJ6/vYVf3AMXeVD59djU3rKqyzv4wFvTNlFJVNnf4eHx9O09saKe9bzgU6D+4pIyLa4vITLUuIjN9AgHlxXe7uO+VXbz8bjfpyYl87PRK/vr8eZTnpEe6ehFnQd9MiX2HBnl8Qzu/eaeN7Qf6SUoQzltYwJXLylhdV2xj501EbOnwce/Lu/jt+jYArlxWxt9dMP+oC8nFCwv65qT5hsd4amMHj77Txlu7DgGwqjqPq5aX8YHGUptNaaJGW+8Q9768k4fe2sfw+AQfXFzKZy9aSE1J/AV/C/rmfZkIKK9s7+bhda38vrmTkfEA8wsz+ciKCq5cWkZlXkakq2jMMR3sH+FHr+zigdd2MzA6wQcWl/D51YtYWBw/wd+Cvjkhu7sHeHhdK4+83UpH3zDZ6clcubSMj55WwdIKmyBjZpfewVF+9Mou7ntlF4NjE1y9rJzPXbyQ6oLYn/hnQd8c0/DYBE9v6uShNXt5Y+chEgTOX1TIX5xWyep6mwpvZr9DA6P84KUdPPDabsYnlBvOqOLWixdSkJUa6apNGwv65k+0dPp58K29PPp2K77hcaryMvjY6ZV8dEUFJdk2lt7EngO+Ye54/l0efGsfaUkJ/PX587jp/HkxucyzBX0DOFf1v9vYwc/f3MPbe3tJSUzgssYSrl9VyZlz822TCxMXdnb185+/b+Gppk5KvGl8+bIarlleHlP//i3ox7ld3QP8/I09PPx2K72DY8wrzOSGVVV8ZEWFrUNv4tba3Yf4xu+2sGFfL4vLs7ntygZOm5Mb6WpNCQv6cWgioDy3ZT8/fWMPL7/bTVKCcFlDCR8/s4qz5uVbp6wxOJO8ntjYzr8+tZVO3zAfWV7OV66opcg7u1OcFvTjSM/AKA+t2cfP3thDW+8QJd40bjijiutWVdq6N8Ycw8DIOHe/sJ0fvrSLlKQEvnDJIj511hySEmfnuj4W9OPA1k4f97+ym8fWtzEyHuDMeXl8+qxqLqkvnrX/cI2Zabu7B7jtiWZeaOmivtTLN69pZEXV7Ev5WNCPUYGA8tzWA9z3yi5e33mQtOQErllewafPnkNtiTfS1TNmVlJVnt7UyT8/sZn9/mE+eeYcvnxZzaxaZsQ2UYkxAyPjPLyulftf3cXug4OUZafxD5fXcv2qSltp0JhTJCJcsbiU8xYV8n9/38KPX9vNs5v3842rGlldXxzp6k2pmdgjtxL4CVACBIB7VPX2Y5W3K/3DdfYN8+PXdvOLN/fgGx5neVUOnzl3Lpc3lFgKx5hp8s7eHr7ySBMt+/1cvayMr324IerXnIqa9I6IlAKlqvq2iHiAdcDVqrr5aOUt6Du2dPj44cs7eXx9OwFVrmgs5TPnzZ2VuUZjZqPR8QB3v7CdO5/fTk5GCt+8upHLG0siXa1jipr0jqp2AB3u734R2QKUA0cN+vFMVXl9x0H+66WdvLSti4yURD5x5hw+c+5cW/DMmBmWkpTA369exKX1JXzp1xv425+t45rl5dx2ZcOs3rd3RnP6IlINLAfenMnPjXYTAacT6b9e3EFTWx8FWal86dJFfOLMOZavNybC6su8/PaWc7jz+e3c+cftvLHzIP/x50s5d2FBpKt2UmZs9I6IZAEvAv+iqo8e8dxNwE0AVVVVp+3Zs2dG6hRpI+MT/ObtNn7w0k52dQ9QnZ/BTefP5yMryklLtkXPjIk2G/b18vlfrWdn1wD/45y5/MMVNVGzQGHU5PTdyiQDTwLPqOp33qtsPOT0B0fH+cWbe/nhyzvZ7xuhsdzL//yzBVzeWEJiDK0FYkwsGh6b4FtPbeEnr++hrtTLHdcti4p1+6Mm6Isz9/8B4JCq/v3xysdy0O8bGuOB13Zz/6u76Bkc46x5+fzdhfM5d0GBLZFgzCzz3Jb9fPnhjQyMjPO1Dzdw/arKiP4/jqagfy7wMtCEM2QT4J9U9amjlY/FoH9oYJQfvbKTn7y2B//IOBfXFnHzRQtsJI4xs9wB/zBf/NUGXn63mw8vLeNb1zRGbEJXNI3eeQWIy8vYLv8IP3x5Jz99fQ/D4xNc0VjCzRcuoKEsO9JVM8ZMgSJPGg/85Sq+/+IOvvPsNja29nLXDStoLI/e/+O2DMM0OOAf5gcv7uTnb+5hdDzAlUvLuPnCBVGR9zPGTI81uw/x2V+8w6HBUb5xVQMfO71qRj8/aq7048kB3zD/5Qb78YBy9bJybrloAXPjYH9OY+Ld6dV5/O7Wc/ncQ+v5h0eaWLu7h29c3Rh1I/Es6E+B4JX9z95wgv1HljvBfk6+BXtj4kl+VioP/I9V3P6Hbdzx/HY2d/j4wSdPoyI3eiZXWnrnFBzsH+EHL+3kJ6/vZmxCuWZ5ObdcuIBqu7I3Ju49v3U/n3toPUkJwl03rODsBdM7mStqRu+8X7Mh6PcNjnHPyzu4/9XdDI9NcPWycj578UJL4xhjDrOre4CbfrKWHV39/NMH6vjMuXOnbVin5fSnQf/IOPe9sosfvrwT//A4H1pSyt+vXsSCoqxIV80YE4XmFmTym5vP4Yu/Ws83f7eFlk4/37ymMaKzeC3on4DhsQl+9sYe7n5hB4cGRrmkvpgvXLKIulLbtMQY896yUpP4/sdP43tunn9X9wDf/8RpFHpSI1IfC/rvYWwiwK/XtnLHc+/S6RvmvIUFfPHSGpZV5kS6asaYWSQhQfjCpTUsKvHwpV9v4Oq7XuX+vzydRREYxm1B/ygCAeV3TR1859lt7OoeYEVVDt/52FLOnj87V9UzxkSHDy0poyovg888sJaP3v0ad39iBectLJzROtjWS2FUlZe2dfHhO1/hsw++Q0piAvd+aiWP/M+zLeAbY6bEkoocHrv5HMpz07nx/jU89NbeGf18u9J3bWzt5dv/vZXXdhykIjed71y7lKuWlduql8aYKVeek86v//YsbvnFO3zl0Sbae4f4/CWLZmTBtrgP+nsODvDvz7Twu40d5GWm8LUP13PDGVVRs0a2MSY2edKSuffTK/nqbzZxx/Pbaesd5tsfXUzyNO99HbdB/2D/CP/v+e387I09JCcmcOtFC/jr8+dFbIU8Y0z8SU5M4NsfXUxZTjrf/cM2uvpHuP/G06c1wxB3QX9odIL7Xt3Ff72wg8GxCa5dWcnnVy+kyJsW6aoZY+KQiPC51QspzUnDNzQ27SnluAn6gYDy6Dtt/OczLXT6hlldV8xXrqhhQZGtfGmMibxrV1bOyOfERdB/bXs33/zdFjZ3+Fhakc3t1y3jjHn5ka6WMcbMuJgO+tsP9POvT23hua0HKM9J547rl/OhxaUk2IgcY0ycmpGgLyKXA7cDicC9qvrt6fy8QwOj3P6Hbfzszb1kJCfylStqufHs6qhb19oYY2batAd9EUkE7gIuAVqBNSLyuKpunurPGh0P8JPXd3P7c+8yODrB9asq+fzqReRnRWaNC2OMiTYzcaW/CtiuqjsBROQh4CpgSoP+3oODfOq+N9l9cJA/W1TIVz9YZ9sTGmPMEWYi6JcD+8IetwJnhBcQkZuAmwCqqk5uX8nSnDRqS7zcdmUDF9QUnWRVjTEmts1E0D9ar+lhO7eo6j3APeBsonIyH5KcmMB/ffK0k3mpMcbEjZlYcK0VCB+AWgG0z8DnGmOMOcK0b5coIknANuBioA1YA9ygqs3HKN8F7DmFjywAuk/h9bNRPJ4zxOd5x+M5Q3ye9/s95zmqetx1mqc9vaOq4yJyC/AMzpDN+44V8N3yp7S4tIisPZF9ImNJPJ4zxOd5x+M5Q3ye93Sd84yM01fVp4CnZuKzjDHGHJttomKMMXEkFoP+PZGuQATE4zlDfJ53PJ4zxOd5T8s5T3tHrjHGmOgRi1f6xhhjjsGCvjHGxJGYCfoicrmItIjIdhH5SqTrM11EpFJE/igiW0SkWUQ+5x7PE5FnReRd98/cSNd1qolIooi8IyJPuo/nisib7jn/UkRSIl3HqSYiOSLysIhsdb/zs2L9uxaRz7v/tjeJyIMikhaL37WI3CciB0RkU9ixo3634rjDjW8bRWTFyX5uTAT9sJU8rwDqgetFpD6ytZo248AXVbUOOBO42T3XrwDPqepC4Dn3caz5HLAl7PG/Ad91z7kH+ExEajW9bgeeVtVaYCnO+cfsdy0i5cCtwEpVbcSZ23Mdsfld/xi4/Ihjx/purwAWuj83Ad8/2Q+NiaBP2EqeqjoKBFfyjDmq2qGqb7u/+3GCQDnO+T7gFnsAuDoyNZweIlIBfBC4130swEXAw26RWDxnL3A+8CMAVR1V1V5i/LvGmT+U7s7mzwA6iMHvWlVfAg4dcfhY3+1VwE/U8QaQIyKlJ/O5sRL0j7aSZ3mE6jJjRKQaWA68CRSragc4DQMQa0uNfg/4X0DAfZwP9KrquPs4Fr/zeUAXcL+b1rpXRDKJ4e9aVduA/wT24gT7PmAdsf9dBx3ru52yGBcrQf+4K3nGGhHJAh4B/l5VfZGuz3QSkQ8BB1R1XfjhoxSNte88CVgBfF9VlwMDxFAq52jcHPZVwFygDMjESW0cKda+6+OZsn/vsRL042olTxFJxgn4P1fVR93D+4O3e+6fByJVv2lwDnCliOzGSd1dhHPln+OmACA2v/NWoFVV33QfP4zTCMTyd70a2KWqXao6BjwKnE3sf9dBx/pupyzGxUrQXwMsdHv4U3A6fh6PcJ2mhZvL/hGwRVW/E/bU48Cn3d8/Dfx2pus2XVT1H1W1QlWrcb7b51X148AfgT93i8XUOQOoaiewT0Rq3EMX4+w4F7PfNU5a50wRyXD/rQfPOaa/6zDH+m4fBz7ljuI5E+gLpoHeN1WNiR/gAzhLOO8A/nek6zON53kuzm3dRmC9+/MBnBz3c8C77p95ka7rNJ3/BcCT7u/zgLeA7cCvgdRI128azncZsNb9vh8DcmP9uwb+GdgKbAJ+CqTG4ncNPIjTbzGGcyX/mWN9tzjpnbvc+NaEM7rppD7XlmEwxpg4EivpHWOMMSfAgr4xxsQRC/rGGBNHLOgbY0wcsaBvjDFxxIK+McbEEQv6xhgTRyzoGxNGRB4TkXXueu43ucc+IyLbROQFEfmhiNzpHi8UkUdEZI37c05ka2/M8dnkLGPCiEieqh4SkXSc5T0uA17FWfPGDzwPbFDVW0TkF8DdqvqKiFQBz6izz4ExUSvp+EWMiSu3isg17u+VwCeBF1X1EICI/BpY5D6/Gqh3logBwCsiHnX2OTAmKlnQN8YlIhfgBPKzVHVQRF4AWoBjXb0nuGWHZqaGxpw6y+kbMykb6HEDfi3OdpQZwJ+JSK67tO9Hw8r/Hrgl+EBEls1obY05CRb0jZn0NJAkIhuBbwBvAG3At3B2J/sDzjK/fW75W4GV7kbVm4G/nfkqG/P+WEeuMcchIlmq2u9e6f8GuE9VfxPpehlzMuxK35jju01E1uOs774LZ117Y2Ylu9I3xpg4Ylf6xhgTRyzoG2NMHLGgb4wxccSCvjHGxBEL+sYYE0f+PyXKXH5mCc12AAAAAElFTkSuQmCC\n",
                        "text/plain": [
                            "<matplotlib.figure.Figure at 0x18b9110c128>"
                        ]
                    },
                    "metadata": {},
                    "output_type": "display_data"
                }
            ],
            "source": [
                "# plot steady state values for consumption, savings and labor\n",
                "plt.subplot(3, 1, 1)\n",
                "plt.plot(chist)\n",
                "plt.title('Consumption', y=.92)\n",
                "plt.xlabel('age')\n",
                "plt.xticks([])\n",
                "\n",
                "plt.subplot(3, 1, 2)\n",
                "plt.plot(nhist)\n",
                "plt.title('Labor', y=.92)\n",
                "plt.xlabel('age')\n",
                "plt.xticks([])\n",
                "\n",
                "plt.subplot(3, 1, 3)\n",
                "plt.plot(khist)\n",
                "plt.title('Capital', y=.92)\n",
                "plt.xlabel('age')\n",
                "\n",
                "plt.show()"
            ]
        }
    ],
    "metadata": {
        "celltoolbar": "Raw Cell Format",
        "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.6.4"
        }
    },
    "nbformat": 4,
    "nbformat_minor": 2
}