Structural Estimation Series III

structural-estimation-series-III.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Structural Estimation Series III のコピー",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyOfkDDULvSlCjLiLbau8jV/",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/gist/justinian336/2677bce43ec4434e7c4da48a4adffda3/structural-estimation-series-iii.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kyCg_ghoyeT5",
        "colab_type": "code",
        "outputId": "d302b999-9ec0-46e1-b1e0-7468cccb7823",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 51
        }
      },
      "source": [
        "import numpy as np\n",
        "from scipy.special import softmax\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "sns.set()"
      ],
      "execution_count": 1,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.6/dist-packages/statsmodels/tools/_testing.py:19: FutureWarning: pandas.util.testing is deprecated. Use the functions in the public API at pandas.testing instead.\n",
            "  import pandas.util.testing as tm\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kM8Qw4838gIL",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Define the DDC Problem:\n",
        "\n",
        "# Number of choices\n",
        "n_choices = 2\n",
        "\n",
        "# Number of states\n",
        "n_states = 5\n",
        "\n",
        "# The discount factor\n",
        "discount_factor = 0.99\n",
        "\n",
        "# The transition matrix:\n",
        "transition_matrix = np.array(\n",
        "        [\n",
        "            [\n",
        "                [1., 0., 0., 0., 0.],\n",
        "                [0.1, 0.9, 0., 0., 0.],\n",
        "                [0., 0.1, 0.9, 0., 0.],\n",
        "                [0., 0., 0.1, 0.9, 0.],\n",
        "                [0., 0., 0., 0.1, 0.9]\n",
        "            ],\n",
        "            [\n",
        "                [1., 0., 0., 0., 0.],\n",
        "                [1., 0., 0., 0., 0.],\n",
        "                [1., 0., 0., 0., 0.],\n",
        "                [1., 0., 0., 0., 0.],\n",
        "                [1., 0., 0., 0., 0.]\n",
        "            ]\n",
        "        ]\n",
        ")\n",
        "# The first element of this array (transition_matrix[0]) represents the transition function for the case when *not replacement* (choice 0) is chosen.\n",
        "# The second element is the stochastic matrix for the case when choice 1 (replacing) is chosen."
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HpvxcjeK84sF",
        "colab_type": "code",
        "outputId": "c22c2ae1-e7aa-4ddd-d4fe-c100ce2ca0f2",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        }
      },
      "source": [
        "# Make sure that the shape of the transition matrix is (2, 5, 5):\n",
        "transition_matrix.shape"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(2, 5, 5)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 3
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qgYv-yyc9rmZ",
        "colab_type": "code",
        "outputId": "91719c06-fd4f-43ad-b8cc-d6f4ba77b6d4",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 51
        }
      },
      "source": [
        "# Make sure that the sum of values across the columns of each matrix (axis 2) equals to 1:\n",
        "transition_matrix.sum(axis=2)"
      ],
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[1., 1., 1., 1., 1.],\n",
              "       [1., 1., 1., 1., 1.]])"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 4
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UHbicMJ7yhdB",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Define the utility function: it returns the utility level for each combination of choices and states\n",
        "# Its output must be a (2, 5, 1) array.\n",
        "def utility_function(theta, n_choices, n_states):\n",
        "    # theta contains the utility function parameters. R is the replacement cost of the \n",
        "    # engine, and mile_cost is the additional mantainance cost\n",
        "    # per mile run.  \n",
        "    [mile_cost, replacement_cost] = theta\n",
        "\n",
        "    # This line produces vectors containing all the possible combinations of states and choices\n",
        "    miles = np.arange(n_states)\n",
        "\n",
        "    # Here we calculate the utility for each possible combination of action and state:\n",
        "    # If action is *not replacing*, the utility is the cost per mile times the number of \n",
        "    # miles since the last replacement\n",
        "    u_no_replacement = mile_cost*miles\n",
        "\n",
        "    # If the replacement is performed, the utility is just the replacement cost, no matter \n",
        "    # the state. We use np.repeat to create an array containing the replacement cost for\n",
        "    # each possible state.\n",
        "    u_replacement = np.repeat(replacement_cost, n_states)\n",
        "\n",
        "    # Finally, the utility is an array \n",
        "    u = np.array([u_no_replacement, u_replacement]).reshape((n_choices, n_states, 1))\n",
        "\n",
        "    # Finally, we reshape the utility function and return. Remember that U is a (n_choices, n_states, 1) tensor.\n",
        "    return u"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5neXmvfD-TIu",
        "colab_type": "code",
        "outputId": "4606d0b8-585e-45a3-c5e9-40db796e136b",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        }
      },
      "source": [
        "# Let's see if the shape is correct...\n",
        "cost_per_mile = -1\n",
        "replacement_cost = -2.5\n",
        "some_theta = [cost_per_mile, replacement_cost]\n",
        "utility = utility_function(some_theta, n_choices, n_states)\n",
        "utility.shape"
      ],
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(2, 5, 1)"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 6
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EBC0cft1-f4p",
        "colab_type": "code",
        "outputId": "e534934d-caa2-4a62-f263-01288183c3bd",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 204
        }
      },
      "source": [
        "# Now let's see its output:\n",
        "utility"
      ],
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[[ 0. ],\n",
              "        [-1. ],\n",
              "        [-2. ],\n",
              "        [-3. ],\n",
              "        [-4. ]],\n",
              "\n",
              "       [[-2.5],\n",
              "        [-2.5],\n",
              "        [-2.5],\n",
              "        [-2.5],\n",
              "        [-2.5]]])"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 7
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lhlTdQQJykpa",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Let's define a function that obtains the discounted value for each combination of choice\n",
        "# and state. It calculates V = U + \\beta*Fv, where v is the **expected discounted value**.\n",
        "# Note that v has dimensions (n_states, 1).\n",
        "# This function is necessary to implement Phi Map and returns a (2, 5, 1) array.\n",
        "def get_discounted_value(utility, discount_factor, transition_matrix, v):\n",
        "    n_choices, n_states, _ = transition_matrix.shape\n",
        "    discounted_value = utility + discount_factor*np.array([transition_matrix[i].dot(v) for i in range(n_choices)])\n",
        "    return discounted_value"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5c_teiIwyoIv",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Phi Map maps the conditional choice probabilities (ccp) into the Value space.\n",
        "# Its output is a (5, 1) array.\n",
        "def phi_map(p, transition_matrix, theta, utility_function, discount_factor, n_choices, n_states):\n",
        "  n_choices, n_states, _ = p.shape\n",
        "  utility = utility_function(theta, n_choices, n_states)\n",
        "  denominator = np.identity(n_states) - discount_factor*((p*transition_matrix).sum(axis=0))\n",
        "  denominator = np.linalg.inv(denominator)\n",
        "  # In the following line I'm using np.nan_to_num to replace any possible ln(0) with zeros.\n",
        "  # Don't worry, it doesn't affect convergence\n",
        "  numerator = (p*(utility + np.euler_gamma - np.nan_to_num(np.log(p), 0))).sum(axis=0)\n",
        "  v = denominator.dot(numerator)\n",
        "  return v"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Q9uYVxH4yqlp",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Finally, Lambda Map is a function that maps from the Value space into the Probability Space.\n",
        "def lambda_map(v, theta, discount_factor, transition_matrix):\n",
        "  utility = utility_function(theta, n_choices, n_states)\n",
        "  V = get_discounted_value(utility, discount_factor, transition_matrix, v)\n",
        "  return softmax(V, axis=0)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dOyIhrj3zGHU",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Let's see the fixed point into action!\n",
        "# Define some random values for V. I'm not fixing the seed to convince you that the fixed\n",
        "# point is reached regardless of the initial values.\n",
        "# You'll see different values, but that's the point (＾▽＾)\n",
        "initial_v = np.random.uniform(size=(n_choices, n_states, 1))\n",
        "\n",
        "# I also initialize the vector of conditional choice probabilities with random values\n",
        "initial_p = np.random.uniform(size=(n_choices, n_states, 1))\n",
        "# Make sure that the probabilities for every choice add up to one in each state.\n",
        "initial_p = (initial_p / initial_p.sum(axis=0))\n",
        "\n",
        "# We will iterate many times using Phi Map and Lambda Map to see how Values and\n",
        "# CCPs converge.\n",
        "repetitions = 5\n",
        "\n",
        "# Let's store the values of the squared deviations of each iteration to visualize the fixed\n",
        "# point magic\n",
        "root_squared_deviations = np.empty(shape=(2, repetitions))\n",
        "\n",
        "p = initial_p\n",
        "v = (p*initial_v).sum(axis=0)\n",
        "\n",
        "some_other_theta = [1,0]\n",
        "for i in range(repetitions):\n",
        "  # P -> V\n",
        "  new_v = phi_map(p, transition_matrix, some_other_theta, utility_function, discount_factor, n_choices, n_states)\n",
        "\n",
        "  # V -> P\n",
        "  new_p = lambda_map(new_v, some_other_theta, discount_factor, transition_matrix)\n",
        "  root_squared_deviation_v = np.sqrt(((new_v - v)**2).sum())\n",
        "  root_squared_deviation_p = np.sqrt(((new_p - p)**2).sum())\n",
        "  \n",
        "  root_squared_deviations[0, i] = root_squared_deviation_v\n",
        "  root_squared_deviations[1, i] = root_squared_deviation_p\n",
        "\n",
        "  p = new_p\n",
        "  v = new_v"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pB4_UNzyNI7T",
        "colab_type": "code",
        "outputId": "913444a5-727f-4a90-b85b-f31fee75ee4b",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 284
        }
      },
      "source": [
        "# Let's visualize now the squared deviations to see if the values of V and P converged:\n",
        "_, ax = plt.subplots(1, 2, sharex=True)\n",
        "\n",
        "for i, title in enumerate(['Squared Deviations of V', 'Squared Deviations of P']):\n",
        "  ax[i].plot(np.arange(repetitions), root_squared_deviations[i])\n",
        "  ax[i].set_title(title)"
      ],
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAELCAYAAAA4HCbKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxcZb348c9M9mWSJmnSSdNlSpcvBVr27QpylU1RrhtbFXvRK4p6QfT6U9y5clWuFxURsCjqRUBAFgG1iOIVsexLS6HAt4U2TZe0SZMuSdOk2X5/nDPtJJkkM8msZ77v16uvzJxz5pxnTp/5zjPPOd/n8Q0ODmKMMcbb/OkugDHGmOSzYG+MMTnAgr0xxuQAC/bGGJMDLNgbY0wOsGBvjDE5wIJ9konIP4vI5jQd+1QR0Um8fpmIfCORZZoMcawSkQ4RuSLd5cl1VrcTJxV1Oz8ZO00GETkF+D5wONAPvA5cqarPp7VgkyQig0AXMAj0AKuAn6nqPZPdt6r+A5AYy3EJ8AlVPSXi9ZdNtgwJ9iXgb6p61PAVIrIMKFXVpcOWHwk8B9Srantqihkfq9vxy6W6DSAijwMnAX1AN/AE8FlVbY71AFnRsheRCuAPwE+AaqAB+E+cCpTqsuQlYbdHqmo5TuX9X+BGEflWEo6T7WYDa0ZZdxvwQREpG7b8o8AfMjjQW902MHbdDvt391wuAKYAP4rnANnSsl8AoKp3uc/3AX8Or3Qr6X8DlwB7gB8ANwIFqtonIo043+yPudtfDcxT1Yvd5/cCpwIlwMvAp1V1jbvuf93jzQZOA94nIq/hfDjfDnQCP1LVG9ztS4CfAu8DmoFfxfomVXUHcLuI7APuEJEbVbVNRCqBHwLnAAPuPr+F8/+3HThFVV91j18LNLnlPQy4Q1VnuOuuAi4F6oBNwNdU9XcishBYBhSISCfQp6pT3Pe+WVW/7r7+UuDLOEFpBXCZqm511w0Cnwb+A6gF7sSpnIMiMg/4BXAU0Av8VVUvjHYORORfgO/hBL1V7v/F6yLyf+75P0VErgeOUdW1EefuaRHZAnwI+LW7rzzgw265MpXVbavbY9btKOeyXUTuJ856nRUte2At0C8it4nIu0Wkatj6S4H3AkcDxwHnxbn/R4D5OBXlJZz/zEgfBr4DBICngN/jfHAagNOBK0XkbHfbbwFz3X9nA/8aZ1kAHsKp7Ce4z/8X5+fbPJz3eBbOB7wHeABYEvHaC4C/q2pLlP2+hfPBr8RpPd4hIvWq+jpwGfC0qpar6pThLxSRd+JU1AuAemAjcPewzd4LHA8sdrcLn5NrcAJYFTADJ5iMICILgLuAK3E+VMuB34tIoaq+E/gHbutmlA/Dr4HIbpwzgAJ3P5nK6rbV7VjqduS+puI0alaOtd1wWRHsVXUPcApO39/PgVYReVhEprmbXABcr6qb3J/r34tz/79U1Q63gl0NHOm2OMIeUtUnVXUAWATUquq3VXW/qq53y3RRRFm+o6rtqroJuGEC77cX2AFUu+/xHJw+3L1uRf9RxPF+E/EYnA/vb0bZ772qulVVB9x+03Uc/NCN5yPAL1X1Jfc8fQU4WURCEdtcq6q7VLUJ+BtOawecFs9sYLqqdqvqilGOcSHwR1X9i3sOrsNpkf5TjGW8HThNRGa4z5cCv3H3lZGsblvdjrGMADeIyC6cL+Nm4AtxvDZrunFwv6EvARCRQ4E7gOtxvvmn4/x0C9sY637dn8nfAc7H+cYdcFdNBXa7jyP3PRuY7p70sDycb2YmU5aIMhW4ZWl3j1cANIscuB7ljzjG34BSETkR52fvUcDvRtnvUpwKEnIXleO8z1hMx2kZAqCqnSLShtMCbHQXb4vYvsvdPzgXn64BnhORncAPVPWXoxzjwPlS1QER2eQeY1yq2iQiTwAXi8iNwPtxuiMymtVtq9sxukJVb41j+yGyJthHUtU33D63T7mLmoGZEZvMGvaSvUBpxPNgxOMP4/RBnoHzH1sJ7AR8EdtEDg26CdigqvNHKV64LOGLLcPLEov34fy0fQ4oxLlYN1VV+4ZvqKr9IvJbnMCwHediZMfw7URkNk4r7XScn7T9IrKKg+9zvOFPt+J8OMP7KwNqgC3jvRlV3YbTHRG+8+QxEXlCVd+McoxFEcfw4ZzLcY8R4TacvtdmnP+nF+N4bdpZ3T7I6nZiZUWwd1s77wHuUdXNIjITpwI8427yW+AKEfkDTuW/atguVgEXicgjwJE4/Z5/ctcFcCpcG86H5rvjFOc5oENEvozzM3Y/sBAoUedWud8CXxGRZ4Ey4PI43mc18G6cC1b/rapt7vI/Az8Q577gTmAOMENV/+6+9DfAg+57+Noouy/DqfSt7j4/BhwRsX47MMPtQ9wf5fV3AXeJyG9wbg38LvCsqjbG8L7Ox/kQbsYJNoMcbGVG+i1wlYicjnNr2edw/m+eGu8YEe4HbsLpt70pjtelhdVtq9vjHSNRsqLPHugATgSeFZG9OB+EV3GujoPzrf4oTl/WSzgXdiJ9A+ei0k6cIBDZ7/drnJ9XW4DXOPghi0pV+3Eu1hwFbMDpf7wVp9WEu/+N7ro/4/Qjj+dl906BN4FPAJ9X1W9GrF+K0wp6zX0P9+FcSAqX6VmcQDAd54JctHK/hnMnx9M4lX8R8GTEJv+H02LbJiI7orz+MZzzeD9OC28uQ/tTx3I8zv9dJ/Aw8Dm3P3j4MRS4GOci1w7gXODcUT6gUanqXreMMxh5MTITWd22up0SPi9OXuJeWNmAe3tamotjTMJY3TYTlS0te2OMMZNgwd4YY3JATN04IvIgzoWTAZyLKJer6io3UeA2nCvXbcBSVV3nvmbUdcYYY1Ir1pb9v6rqkap6NE4yQPg+0mXATaq6AOfOh1siXjPWOmOMMSkU9wVaN3nhCpzMt7VAjXtfax5OC34+zv2tUdepamsMhynCucrdjDMKoDGJlIdzx8fzpHbAMavXJpnC9boZ566pIRfwY77PXkRuxRm3wge8CzchwL1dK5wAsdVd7htjXSzB/ngOZu0Zkyyn4gx6lSpWr02qzOFg9i8QR7BX1U8AiMhHgf/BuS81WWIeo9mYSUh1PWsG2LlzLwMDI39R19SU09bWmeIiZTY7J9FFOy9+v4+qqjJwGjEjJpWJO4NWVW8XkZ+5O2sQkbyIrprw2Bm+MdbFoh+gra0z6oeitjZAa+uIrOmcZudkpNHOid/vo6amHFLfldIPMDAwGLVeh9eZoeycRDfGednMsC4ciOECrYiUuync4efn4gxi1IKTqh0egnQJsFJVW93R66Kui/F9GGOMSaBYWvZlwL3u4ED9OIH+XHUG7r8MuE1EvomT6hw5lvhY64wxxqTQuMFeVbfjzH0Ybd0bOON6xLXOGGNMalkGrTHG5AAL9sYYkwMs2BtjTA7IqmDf29fPl5c9xZr1bekuijEJ9fwbLXzh+r8z4MEhx01myKpg7/P52NW5n2detZwr4y3d+/tYt2kX29u70l0U41FZFezz8/yEggHeaGxPd1GMSag5wQoANm6zxDiTHFkV7AHmNlTy5ubd9PZFm+bRmOxUP7WUwoI8Gi3YmyTJvmA/vZK+/gGattuHwnhHnt/PIdMraGzek+6iGI/KvmDf4PzcfWvL7jSXxJjEmjdzChtboo8HZcxkZV2wn1JeRF1VCW9utRaQ8ZZ5M6bQs7+fbXaR1iRB1gV7gENnV1vL3njOvJlTAGjcZg0Zk3hZGewlVMXOjh7a93SnuyjGQ0TkOhHZICKDInLEONuKiHSJyHWJOv6MugCFBX4am+16lEm8rAz2h86uBmC9deWYxHoQeDvOlG6jcudnuMXdPmHy/D5mTwvQaDcfmCTIymA/Z3olBfl+3rSuHJNAqrpCVWOZYOcq4A848ywn1OxggKbtHfQP2K3FJrHinqkqExTkO8lVb221YG9SS0SOBM4G3sEEp+Z0Z8mKavGCOh57YTM9A04r3zgzjpmR4j0vWRnswUmueuyFTfT2DVCQn5U/UEyWEZEC4GfAx9zpNie0n7Gm26wpKwBg5evbKM33TaK03mDTbUYX7bxETLcZVdZGSSe5atCSq0wq1QNzgeUi0ghcCVzqzsmcENOqSykqzLOLtCbhsrhlfzC5am5DZZpLY3KBqjYBU8PPReRqoFxVv5ioY/h97kVau/3SJFjWtuynlBdRU1FsyVUmYUTkBhHZDMwAHhORNe7y5SJyXKrKEQoGaGrptIu0JqGytmUPTut+3Wa7SGsSQ1WvAK6IsvycUba/OhnlCNUH6H1+gK07uphZN3ofrDHxyNqWPTgXaS25ynhNyB3u2AZFM4mU1cF+nttX/5Z15RgPqasqoaTIhjs2iZXVwX5mXTkF+X4bJ8d4ysGLtBbsTeJkdbAPz1xlyVXGa0LBCja1dNLXbxdpTWJkdbAHp99+47YOm7nKeEqoPkBf/wBbWvemuyjGI8a9G0dEaoDbcZJJ9gPrgE+paquIDAKvAOFI+1FVfcV93bnA/7jHeBEn6zDhA3XPnV7Jn/qbaNreYffbG88IBZ1U+MZte5gdtOECzOTF0rIfBL6vqqKqi4C3gGsj1v+Tqh7l/gsH+nLg58C5qjoP6AASlngSyWauMl5UO6WE0qJ8m4DcJMy4wV5V21X18YhFzwCzx3nZu4EXVHWd+3wZcOGESjgOS64yXuTz+ZgdDLDBgr1JkLj67EXED3waeDhi8eMiskpEviciRe6yWQwdE7wJmDmpko5hbkOFteyN54TqA2xu6bTrUSYh4s2g/QnQCdzoPp+lqptEpAKnX/8bwNcTVbixRnCLHN7zSKnjuddb8BXkM3VKSaIOn3VsKNiRsvmczAlW0D8wyJYdnQcSrYyZqJiDvTv92nycfvgBgPBED6q6R0RuBb7gbt6EM9532CwglkkhhhhrKNjI4T2DlcUAPPfKVo4/tC7ew3iCDQU70mjnZLyhYDNF+MJsY3OHBXszaTF144jId4Fjgferao+7rEpEStzH+cB5wCr3JX8CjheR+e7zy4DfJrLgkSy5ynjR1MpiyorzbQRMkxDjBnsRORz4CjAdeMrtn/8dcCjwrIi8DKwGenFn7lHVDuCTwB9E5E2gEkjYxMzDWXKV8SKfz0eovsLGtjcJMW43jqquAUabMmfxGK97CHhoguWKm81cZbwoFAzwp2eb6O3rpyA/L93FMVnMM1HRZq4yXhQKBugfGGRTi2XSmsnxTLCfZ8lVxoMODHds/fZmkjwT7CvLi5haaclVxluqK4oIlBbYCJhm0rJ6pqrhDpluM1eZiXNvL/4QEAIWqeqrUbb5BnAR0I9zU8JXVfXRZJUpnElrF2nNZHmmZQ82c5WZtAeBtzM0+3u454DjVXUx8HHgnvAtyMkSClawdcdeenr7k3kY43GeatlHzlxVXVGc5tKYbKOqKwBEZKxtIlvxq3HuVKsBNierXHOCAQYGB9nU0nmgjhsTL0+17C25yqTYUuAtVU1aoIeDmbQ2AqaZDE+17A8kV1mwN0kmIqcB1wBnxvvaWMd8Cps6tZwpgSKad+7L6rF+JioX33Ms4j0vngr2YMlVJvlE5GTgDuB9qqrxvj7WMZ8izaorRxvbc278IxvzKbpo52W8MZ88Fw0tucokk4gcD9wDnKeqL6XquKFggK1te+nZbxdpzcR4LthbcpWZKBG5QUQ2AzOAx0Rkjbt8uYgc5252M1AC3OKOE7VKRBYlu2yhYAWDg9DUYo0YMzGe68aJTK46K92FMVlFVa8Aroiy/JyIx8entFCuyOGO58+Yko4imCznuZY9OP321rI3XlIVKGJKeaENm2AmzJPB/pDpFZZcZTwnFKywYRPMhHky2EcmVxnjFaFggG1tXezr6Ut3UUwW8mSwt+Qq40Wh+gCDYHeamQnxZLC35CrjRbMPDHdswd7Ez5PBHpyLtBu3d9DbN5DuohiTEJVlhVQFimzYBDMh3g32llxlPCgUDLDBgr2ZAM8Ge0uuMl4Uqq9ge3sXXd12kdbEx7PB3mauMl4UcpOr7BeriZdngz1YcpXxngOZtNaVY+Lk7WBvyVXGYypKC6mpKLZMWhM3bwd7S64yHhSqtzlpTfw8Hewtucp4USgYoGXXPvZ296a7KCaLeDrYW3KV8aKQm1xl99ubeIw7xLGI1AC3A3OB/cA64FOq2ioiJwG34Izv3QhcrKot7utGXZdKNnOV8ZrIi7SHharTXBqTLWKJfoPA91VVVHUR8BZwrYj4caZm+6yqLgCeAK4FGGtdqoWTqzbarWrGI8pLCphaWWx35Ji4jBvsVbVdVR+PWPQMMBs4FuhW1RXu8mXABe7jsdalVDi5ar115RgPCdVX0NhsNx6Y2MXVr+G22D8NPAzMAjaG16nqDsAvItXjrEspS64yXjQnGGDH7m4699lFWhObeKcl/AnQCdwIfCDxxRlqrJnSa2sDMe/nsENqeG19W1yvyUZef38TEc85EZHrgA8BIWCRqr4aZZs84AbgXThdnNeq6q0JKWwcQgf67fdwxJyaVB/eZKGYg737QZgPnKuqAyLShNOdE14/FRhQ1fax1sVTuLa2TgYGBkcsr60N0Noae3/ljJpSnli5BX2rleqK4niKkDXiPSe5YLRz4vf7RmtIPAj8GPjHGLv9CDAP57NQA6wUkcdUtXHSBY5D+CLtxm0dFuxNTGLqxhGR7+L0w79fVXvcxS8CJSJyivv8MuDeGNalnCVXmVio6gpV3TTOZhcCP1fVAVVtxfmCOD/5pRuqtLiAuqoSS64yMRs32IvI4cBXgOnAUyKySkR+p6oDwEeBn4rIOuA04CqAsdalgyVXmQQacj0KaAJmpqMgoWDAhk0wMRu3G0dV1wC+UdY9BSyKd12qWXKVySSJuhZ1+Nxannu9hcKSQirLixJRtIxk16Kii/e8xHuBNmtZcpVJkPD1qOfd58Nb+uNK1LWo2kAhAC+uaWbRId7st7drUdFFOy9jXIty1ie7UJnCkqtMgtwLXCoifhGpBd4P3JeOghzIpLX77U0McibY28xVZjwicoOIbAZmAI+JyBp3+XIROc7d7HZgPc6wIc8A31bVDekob0lRPsHqUsukNTHJmW6ccHKV3ZFjRqOqVwBXRFl+TsTjfpzEwowQCgbQTbvSXQyTBXKmZQ82c5XxnlAwwM6OHnZ39oy/sclpuRXsbeYq4zGheqd70rpyzHhyK9hbcpXxmFnTyvFhwd6ML6eCvSVXGa8pLswnWFNqE5mYceVUsLfkKuNFoWAFGyyT1owjp4I9OF05G7d30Ns3kO6iGJMQofoAuzv3s7PDLtKa0eVesLfkKuMxoYgRMI0ZTc4Fe0uuMl4zqy6Az4cNimbGlHPB/kBylQV74xFFhXlMn1pmd+SYMeVcsAc3ucpuvzQeEpoWoLF5D4ODIwdYMwZyNdhbcpXxmFB9BXu6eu0irRlVbgZ7S64yHnNwTlrryjHR5WSwt+Qq4zUz68rx+3x2kdaMKieDvSVXGa8pLLCLtGZsORnswZKrjPeE6gM0NnfYRVoTVe4Ge0uuMh4zJxigc18vbXbjgYkiZ4O9JVcZr5kddIc7brYGjBkpZ4O9JVcZr5lZV0ae32e/Vk1UOTMtYTRzGypZa1O6GZeILABuA2qANmCpqq4btk0d8CtgJlAA/A24QlX7UlzcEQry82ioLbMJyE1UOduyB0uuMiMsA25S1QXATcAtUbb5KvC6qi4GFgPHAh9MXRHHFgoGaNxmF2nNSLkd7C25yrjcFvsxwF3uoruAY0Skdtimg0BARPxAEVAIbElZQccRClawt7uPHbutAWOGiqkbR0SuAz4EhIBFqvqqu7wR6Hb/AXxZVR91152E0zIqARqBi1W1JXFFn7yZdeUUuslVxx9al+7imPSaCWxR1X4AVe0Xka3u8taI7a4B7geagTLgRlV9Mp4D1dSUj7qutjYQZ7GHOnphkF8/qrR39XLYfG/U6cmeE6+K97zE2mf/IPBj4B9R1p0XDv5hbqvnDuASVV0hIl8HrgU+HlfpksySq8wEnA+sBk4HAsAjInKeqt4X6w7a2joZGBjZzVJbG6C1dXIXV0vzfeTn+VitLcj0ikntKxMk4px4UbTz4vf7xmxIxNSNo6orVHVTHGU5FuhW1RXu82XABXG8PmUOseQq49gENIhIHoD7d7q7PNLlwJ2qOqCqu4GHgHektKRjKMj301Bbbpm0ZoRE9NnfKSKrReRmEZniLpsFbAxvoKo7AL+IVCfgeAllyVUGwO1iXAUscRctAVaqauuwTTcA7wIQkULgDOBVMsgcu0hropjsrZenquomESkCrgduBC6efLEcyezbDDuxqICbfvcK23d3c/JRMxKyz3Sxvs2R4jwnlwG3icg3gZ3AUgARWQ58U1VfAK4ElonIK0Aezq2XP09ooScpVF/B46u20rJrH9OqStNdHJMhJhXsw107qtojIjcDD7urmoDZ4e1EZCowoKrt8ew/mX2bkaZWFvOytvC2w6YlbJ+pZn2bI412Tkbr21TVN4AToyw/J+LxW8CZiS1pYs2e5g533Nxhwd4cMOFuHBEpE5FK97EPuAjnZzDAi0CJiJziPr8MuHcyBU0mm7nKeElDbRn5eX6bgNwMEVOwF5EbRGQzMAN4TETWANOAx0VkNU6f5QLgMwCqOgB8FPipiKwDTgOuSkL5E8KSq4yX5Of5mVlXbmPbmyFi6sZR1SuAK6KsOnqM1zwFLJpguVIqMrmquqI4zaUxZvJC9QGefnUbA4OD+H2+dBfHZICczqANi0yuMsYLQtMCdO/vZ3t7V7qLYjKEBXssucp4T6jeSaiyfnsTZsHeZclVxkumTy2lIN9vyVXmAAv2LkuuMl6S5/czq67chjs2B1iwd9nMVcZrQsEKNm6Pnqtico8Fe5fNXGW8JlQfoKe3n212kdZgwX4IS64yXhIKupm0dr+9wYL9EJZcZbykvqaMwgK/TUBuAAv2Q4STq960rhzjAX6/j1nTAjTaTQcGC/ZDhJOr1ltXjvGIUDBA0/YO+gfsluJcZ8E+giVXGa8JBQPs7x2guc0u0uY6C/bDzLXkKuMhoaBzS7H12xsL9sMcYslVxkOC1aUUFebZsAnGgv1wllxlvMTv9zF7WsBuvzQW7Iez5CrjNaFggKaWTrtIm+MmOwetJ81tqGTtpl3pLoZJMRFZANwG1ABtwFJVXRdluwuAbwA+YBA4Q1W3p7Ks8QgFA/T2DbB1Rxcz60af19l4m7Xso7Dkqpy1DLhJVRcANwG3DN9ARI4DrgbOVNUjgFOAjP4ZGB7u2AZFy20W7KOw5KrcIyJ1wDHAXe6iu4BjRKR22KafB65T1W0AqrpbVTO6VVBXVUJxYZ4Nd5zjrBsnioMzV+3hhIXT0l0ckxozgS2q2g+gqv0istVd3hqx3WHABhF5AigHHgC+o6oxDy1ZUzN6V0ptbWACRR/f/JlVbN6xN2n7T6ZsLHMqxHteLNhHcSC5aqu17M0IecBi4EygEPgT0AT8OtYdtLVFH3a4tjZAa2tyWt8NNaU89uJmmrftJj8ve37QJ/OcZLNo58Xv943ZkMie//UUm9tQSZMlV+WSTUCDiOQBuH+nu8sjNQH3qWqPqnYADwEnpLSkExCqD9DXP8CW1r3pLopJEwv2o5jbYMlVuURVW4BVwBJ30RJgpaq2Dtv0N8BZIuITkQLgdODl1JV0YmbbcMc5z4L9KOZOt+SqHHQZcLmIrAUud58jIsvdu3AA7gZagNdwvhzWAL9IQ1njUjelhJKifMukzWHWZz8KS67KPar6BnBilOXnRDweAL7g/ssaPp+PUDDABgv2Octa9mOwmauMl4SCATa3dNp1qBxlwX4MllxlvCRUX0H/wCCbWzvTXRSTBuN244jIdcCHgBCwSFVfdZePmloea9p5potMrjqhojjNpTFmcsJz0m7c1sEcN6vW5I5YWvYPAm8HNg5bPlZq+bhp59kgMrnKmGw3tbKYsuJ8uyMnR40b7FV1haoOudd4rNTyONLOM54lVxkvCV+ktYlMctNE78YZK7XcN8a64fcsjykdaeXDLZpfy0NPvEXllFIKC/JScsyJsrTykeycDBWqr+BPzzbR29dPQX5m12eTWBl962U60sqHq68qoa9/kBfXNDPP7cPPRJZWPtJo52S8tHIvCwUD9A8MsqllL4dMt377XDLRu3HGSi2PNe08K1hylfESy6TNXRMK9mOllseRdp4VLLnKeElNRTHlJQXWb5+Dxg32InKDiGwGZgCPicgad1XU1PIY1mUdS64yXuHz+QjVB2xs+xw0bp+9ql4BXBFledTU8vHWZaO50yt49rXttO/pptrutzdZLhSsYPmGjfT09lOU4TcdmMSxDNoY2MxVxktCwQADg4NsarFM2lxiwT4GllxlvCQyk9bkDgv2MbDkKuMlVYEiKsoKbQLyHGPBPkZzGyrZuK2D3r7+dBfFmEk5kElrLfucYsE+RnMbKukfGGTjNuvnNNkvFAywtW0vPfut8ZIrLNjH6EBylXXlGA8IBSsYHISmFmvd54qMHi4hk1hylffFMzS3iAiwErhZVb+YulImxoFM2uYO5s+YkubSmFSwln0cLLnK82IamtsdAuQWnOG/s1JVoIjK8kIbNiGHWLCPg81c5V1xDs19FfAHYG2KipcUc4IVdpE2h1iwj4MlV3naiGG7gfDQ3AeIyJHA2cCPUl7CBAsFA2xr62JfT1+6i2JSwPrs4xCZXHXCwmnpLo5JMREpAH4GfMydp2FC+8mEeRoAFksdD67YwJ6efmbNqErZceNlcxJEF+95sWAfB0uu8rQDQ3O7gTza0Nz1wFxguRvopwA+EalQ1U/GeqBMmKcBoKrE+fivemM70yqKUnbceNg8DdFFOy/jzdNgwT5Ocxsq+fPzm2ymH49R1RYRCQ/NfQdRhuZW1SZgavi5iFwNlGfj3Tjg3GFWFSiyYRNyhPXZx8mSqzwt6tDcIrJcRI5La8mSJBQMsMGCfU6wln2cIi/SzpuRudMUmviNNjS3qp4zyvZXJ7tMyRYKBli5bgdd3X2UFls48DJr2cepsqyQqZXFrLd+e+MBoXonM3zjdmvde50F+wmw5CrjFbNtuMHc2osAABAkSURBVOOcYcF+Aiy5ynhFRWkhNRVFlkmbAyzYT4AlVxkvCQUrbALyHGDBfgJs5irjJaH6AC279rG3uzfdRTFJZMF+Aiy5yniJ9dvnBgv2E7QwVM36rXtYsbo53UUxZlJCQeeOHBsUzdss2E/QOSfN5vA51fzqkdd54Y2WdBfHmAkrLylgamWxzUnrcRbsJ6gg38+/f2ARcxsqueXhNbyyvi3dRTJmwkL1Ntyx11mwn4SiwjyuPO9IZtSWc9MDr6BNO9NdJGMmJBQMsGN3N5377CKtV006P1pEGoFu9x/Al1X1URE5CWc2nxKgEbhYVT3X31FanM/nLzyS/77zJX5832r+35KjmeNmJRqTLULhaQq37eGIOTVpLo1JhkS17M9T1aPcf4+KiB9n5MDPulO8PQFcm6BjZZyK0kK+eNHRlJcU8MN7VrGl1QZJM9nF7sjxvmR14xwLdKvqCvf5MuCCJB0rI1QFivjikqPJz/dz3T2raNnZle4iGROzsuIC6qaUWHKVhyUq2N8pIqtF5GYRmQLMAjaGV6rqDsAvItUJOl5GqptSwhcvPIr+/kGuu3uVDadgskqoPmDDJnhYIsY0PVVVN4lIEXA9cCPwuwTsN2Omb4tHbW2Aaz5Vwld/+iTX37eaaz97CpXlqZkFKFPPSTrZOYldKFjBc6+3sKdrPxWlhekujkmwSQd7Vd3k/u0RkZuBh4EfA7PD24jIVGBAVdvj2XemTN8Wr8riPD533mJ+eM8qvnrzCr605GhKiwuSesxMPyfpMNo5GW/6tlwViui3X3SIXaT1mkl144hImYhUuo99wEXAKuBFoERETnE3vQy4dzLHyjYLZk7hsx9cxJbWvVx/72p69venu0jGjGnWNPeOHEuu8qTJ9tlPAx4XkdXAq8AC4DOqOgB8FPipiKwDTgOumuSxss6iQ2r41L8czltbd3PjA6vp7RtId5GMGVVpcT7TqkstucqjJtWNo6rrgaNHWfcUsGgy+/eC4w6t42P7F/LL5a+z7KFX+cwHjiDPb7lsmUhEFgC3ATVAG7BUVdcN2+YbOL9g+4Fe4Kuq+miqy5osc4IBdNOudBfDJIFFnRQ4ZXE9Hz5jPivX7eCXf3yDgcGR1yFMRlgG3OTmhtyEkxQ43HPA8aq6GPg4cI+IlKSwjEk1OxhgZ0cPuzt70l0Uk2AW7FPkjONm8oG3H8LTa7Zx51/WMmgBP6OISB1wDHCXu+gu4BgRqY3cTlUfVdVwEsVqwIfzS8ATDmbSWleO11iwT6H3njybd584i7+9tIX7/74+3cUxQ80EtqhqP4D7d6u7fDRLgbdUdXMKypcSs6YF8AHPvrbdrjF5TCLuszcx8vl8nPfPc9m3v5/lz2ykpCiP95wcSnexzASIyGnANcCZ8b420/NH3vO2OfzhyQ1sadvL5ecfzcI56c2FzIRzkoniPS8W7FPM5/Nx8VkL6N7fx/1/X09xYT6nHzsj3cUysAloEJE8Ve0XkTxgurt8CBE5GWfsp/epqsZ7oEzPH/ngqXOYNz3Arx9VvnzjP3jHMQ186LS5lBSlPlxkyjnJNNHOy3j5I9aNkwZ+n4+Pn7OQo+dP5c6/rOXJV2y2q3RzR2RdBSxxFy0BVqpqa+R2InI8cA/O4H8vpbaUqbN47lSu+bcTOf3YGfztpS18/dZnWbVuR7qLZSbBgn2a5Of5uex9h7NwdhW/XP46L3pv9OdsdBlwuYisBS53nyMiy0XkOHebm3GG7b5FRFa5/zx5i3FJUT4fPnMBX/3osZQW5XPD/av56YOvsnvv/nQXzUyAL0PvCgkBGzL9524idO/v4wf3rKKxuYPPnbeYIyaYpu6lc5IoMQyXMAdnroVUCZGl9bqvf4BHntnI759qpKggjwveOY9TFtXj8/mSetxMPifpNE43TtR6bS37NCsuzOfz5x9Jw9QybnzgFdZaQovJQPl5fs592xz+8+MnMH1qGb9a/gbX3W1DeWcTC/YZoLS4gC9ceBTVFcX8+L6XbQIJk7Hqa8r48keOYenZQuO2PXzjF8/xyDMb6R+w2zQznQX7DFFRVsgXLzqK0qICfnDPKrbs2JvuIhkTld/n45+PbuC/PnESR8yp5t7H3+Ka216wRkqGs2CfQaorivnikqPI8/v4wd0radm1L91FMmZUVYEi/v2Di/jM+49gd+d+rrntBX77tzfp6bURXjORBfsMM62qlP+46Ch6+wa47q6V7OywMUpM5vL5fBx3aB3/demJnLK4nj8928Q3f/EsrzXGNXWFSQEL9hloRm05X7jwKDr29XLd3Svp6LJb3UxmKysu4JJ3H8qXlhyN3+fjurtX8Ys/vkbnvt50F824LNhnqDn1FVx53mJ27O7mh/e8TFd3X7qLZMy4Dp1dxX9+/ATec/Jsnn51O1//+TM89/p2G/gvA1iwz2Ayq4rPfuAINrd28uP7Xra+UJMVCgvy+NBpc/nmJcdRXVHMsofWcMN9q2nf053uouU0C/YZbvHcqXzyXw7nzS27ufGBV2wkQpM1Zk0L8LWlx3LhO+fxetNOvnbrs/z1xc02n0OaWLDPAscfWscl7zqUNRva+dnDa+yeZpM18vx+zj5hFtf824nMa6jkzr+s5Xt3vMiW1s50Fy3nWLDPEqceOZ2LTp/Pi2tb+dVym+3KZJfaKSV84YIjufS9h7G9fR9X/+p5HvzHevulmkI2xHEWOev4mXT39PHgig2UFObz4TPnJ31sEmMSxefzcfIRQQ4/pJq7/7qOh59s5Pk3WvjYuxcyb0ZluovnedayzzLnvi3E2SfM5K8vbeaBJ2y2K5N9KkoL+eS5h3Pl+Ueyv7ef793xInf8WdnXY3ecJZO17LOMz+fjgnfMY19PP398eiMlRfmcc9LsdBfLmLgtnlvDNZ84kQeeWM9fX9jMynU7uPisBRw9v3b8F5u4Wcs+C/l8PpaeLZywsI77Hn+L/3vJM1OgmhxTXJjPh89YwFeXHktpcT4/uf8Vbn7wVXZ3WuZ4olnLPkv5/T4+8d7D2N87wB1/XsuG7Z3UlBdSX1NGsLqUYHUpRYV56S6mMTGZO72Sb11yPI8828Tvn9zAaxvaufCd8/jA6QvSXTTPsGCfxfLz/Hz6/Ydz+5/Xsm7TLp5q7yLyHp3qiiLqq0sJul8A9TXOl0BVoMgu7JqMk5/n59x/CnGc1HLbI2/wq0fe4JnXWwhNK6e6opiaimKqK4qoqSymtCjf6nCcLNhnuYL8PD5+zkJqawNs2bqLlp37aG7vYlvbXprbu2hu62LFK8307D+YfVtUmEew6mDwD4b/VpdSWGC/Bkx61deU8aWPHMMTq7by15c285cXdtLXP/RW46LCvIPB3/0iiHw+JVBEfp71UkdKarAXkQXAbUAN0AYsVdV1yTxmLissyGNGXTkz6obOMD84OMiuzv0HvgC2tXXR3N7Fus27eea17Qe28+EMsxz+Ejj4ZVDGlPJCz7ekYqmvIpIH3AC8CxgErlXVW1NdVq8Lj5l//lmHsr1lDx1dvbTv6aZtdzfte7rZsaeb9j09tO3pZuO2Djq6hg645gOmBIoOBP/q4V8IOfjrINkt+2XATap6h4hcDNwCvDPJxzTD+Hw+qgJFVAWKWBiqHrKup7ef7e1dbIv4EtjW1sXazVvZ33sw4aW4MG/IF0D42sC06hIK8j3zayCW+voRYB4wH+dLYaWIPKaqjSktaQ7x+3xUlhVSWVbInPqKqNvs7+2nvcMJ/u27u52/7pdB47YOXlrbOu6vA+cL4eAvBa/9OkhasBeROuAY4Ex30V3AjSJSq6qtyTquiU9RQR6zpgWYNS0wZPng4CA7O3oOBH/ni2AvumkXT68Z+mugprL4QPCvnVJMnt9tLfl8+CI3PPhnSIsqcpvwK4Y3uHwR64bva7TXva24MMazEFd9vRD4uaoOAK0i8iBwPvA/MR/MJFxhQd6BrshoBgYHR/w6aNvT4/7tpmlbB3vG+XVQWlyAzwd+fE599Dl//b6hz8N1Ndo653nEMjjw2D9kH1Fej/O3pCiPd9SUR3ubY0pmy34msEVV+wFUtV9EtrrLYwr2NWO8odrawKjrclWiz0ldHcjckcu7e/rY0trJ5pbOg39bOvn7y1vZn0Ejc27e0cWnPrg41s1jra+zgI0Rz5vcbWJm9To+iTon08ZZ39Pbz45d+2jd2UXrzn207trn/u1iU+teurp7GRx0GkID7t8DjwdGLkum2Q1VzJs5Ja7XZPQF2ra2TgainLXa2gCtrTbfZaRUn5OKojwOm1nJYTMPprkPDA7Suc/5QDA4eODOoOHD+ESObR5+OMggHHg89G/kvkZsE3VfzpMjZFrUc+L3+8YMuMlm9Tp2qT4nhUBDVQkNVSWT3tegW2+dL4DhXxIM+VIYiFh2cN3QZeHxsAry/cybOWXEeRmvXicz2G8CGkQkz20l5QHT3eXGg/w+HxWlsXedJFuc/a2x1tcmYDbwvPt8eEvfGOBg18uIPsk0SdrVB1VtAVYBS9xFS4CV1l9vMlEc9fVe4FIR8YtILfB+4L7UldSYiUn2pebLgMtFZC1wufvcmEwVtb6KyHIROc7d5nZgPbAOeAb4tqpuSEdhjYlHUvvsVfUN4MRkHsOYRBmtvqrqORGP+4FPp7JcxiSCd24iNcYYMyoL9sYYkwMs2BtjTA7I1Pvs88C5b3Q0Y63LVXZORop2TiKWpXqcB6vXE2DnJLrh5yXi+QxgMzBk6i/f4PCMl8xwCvCPdBfCeN6pwIoUHs/qtUmVOUBj5IJMDfZFwPFAM5A5+ffGK/KAepzEqFROiWT12iRTuF434yT6ZUXL3hhjTALZBVpjjMkBFuyNMSYHWLA3xpgcYMHeGGNygAV7Y4zJARbsjTEmB1iwN8aYHJCpwyVEJSILgNuAGqANWKqq69JbqvQRkRqc8dXnAvtxxlj/lE0Q4xCRbwFXA4tU9dU0F2dUVq9Hsro9tonU7Wxr2S8DblLVBcBNwC1pLk+6DQLfV1VR1UXAW8C1aS5TRhCRY4CTyI4pA61ej2R1exQTrdtZE+xFpA44BrjLXXQXcIw7NVxOUtV2VX08YtEzOPOj5jQRKcIJmhk/yYjV6+isbkc3mbqdNcEemAlscWcKCs8YtNVdnvNExI9TAR5Od1kywLeBO1S1Md0FiYHV63FY3R5iwnU7m4K9GdtPgE7gxnQXJJ1E5GTgOODmdJfFJIzVbSZft7Mp2G8CGkQkD8D9O91dntNE5DpgPnChqg6kuzxpdhqwENggIo04Y3s/KiJnpbNQY7B6PQar20NMqm5n1aiXIvI4cKuq3iEiFwP/pqrvSHOx0kpEvgucDLxHVbvSXZ5M434o3pvhd+M8jtXrEaxujy3eup1Vt14ClwG3icg3gZ3A0jSXJ61E5HDgK8Ba4CkRAdigqh9Ia8FMvKxeD2N1O/GyqmVvjDFmYrKpz94YY8wEWbA3xpgcYMHeGGNygAV7Y4zJARbsjTEmB1iwN8aYHGDB3hhjcoAFe2OMyQH/H5ZbE/a9kiBIAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 2 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "y_cQZbT3zPoq",
        "colab_type": "code",
        "outputId": "df134657-c4cb-40b2-a8fd-86eb0345bf3e",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 102
        }
      },
      "source": [
        "# It looks like it converged quite fast!\n",
        "# Now let's see if the expected discounted value has converged to a fixed point.\n",
        "# This is v:\n",
        "v"
      ],
      "execution_count": 13,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[127.03628455],\n",
              "       [130.07404485],\n",
              "       [141.78482845],\n",
              "       [161.59553837],\n",
              "       [188.76306388]])"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 13
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wAhtUnyI2LvX",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# And this is the expected discounted value calculated from Equation (9)\n",
        "fv = np.array([transition_matrix[i].dot(v) for i in range(n_choices)])\n",
        "expected_v_calculated = (p*(utility_function(some_other_theta, n_choices, n_states) + np.euler_gamma - np.log(p) + discount_factor*fv)).sum(axis=0)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Z6n4513kK-DD",
        "colab_type": "code",
        "outputId": "f0405c39-f419-48b4-86d0-7a7a85c0af9f",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 102
        }
      },
      "source": [
        "# Let's see if it is the same as v:\n",
        "expected_v_calculated"
      ],
      "execution_count": 15,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([[127.03628455],\n",
              "       [130.07404485],\n",
              "       [141.78482845],\n",
              "       [161.59553837],\n",
              "       [188.76306388]])"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 15
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JQGGhxOtLV_i",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# It looks like it converged!"
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}